{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import batman\n",
    "\n",
    "from lmfit import Parameters, minimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def transit_model(x, epo, per, rprs, ars, b, u1, u2):\n",
    "    '''\n",
    "    Calculates transit model using batman\n",
    "    INPUT:\n",
    "        x - array of times to calculate model on (days)\n",
    "        epo - epoch (days)\n",
    "        per - orbital period (days)\n",
    "        rprs - planet radius divided by stellar radius\n",
    "        ars - distance from planet to star divided by stellar radius\n",
    "        b - impact parameter\n",
    "        u1 - first quadratic limb darkening parameter\n",
    "        u2 - second quadratic limb darkening parameter\n",
    "    OUTPUT: array of model fluxes\n",
    "    '''\n",
    "    bparams = batman.TransitParams()\n",
    "    bparams.t0 = epo\n",
    "    bparams.per = per\n",
    "    bparams.rp = rprs\n",
    "    bparams.a = ars\n",
    "    bparams.inc = np.arccos(b/ars)*180./np.pi\n",
    "    bparams.ecc = 0.\n",
    "    bparams.w = 90.\n",
    "    bparams.u = [u1, u2]\n",
    "    bparams.limb_dark = \"quadratic\"\n",
    "    m = batman.TransitModel(bparams, x)\n",
    "    model = m.light_curve(bparams)\n",
    "    return model\n",
    "\n",
    "def residual(params, x, y, yerr):\n",
    "    '''\n",
    "    Calculates residual between observations and model\n",
    "    INPUT:\n",
    "        params - collection of lmfit parameters\n",
    "        x - array of times (days)\n",
    "        y - array of observed fluxes\n",
    "        yerr - array of uncertainties in observed fluxes\n",
    "    OUTPUT: difference between observed flux and modeled flux\n",
    "    '''\n",
    "    epo = params['epo']\n",
    "    per = params['per']\n",
    "    rprs = params['rprs']\n",
    "    b = params['b']\n",
    "    ars = params['ars']\n",
    "    u1 = params['u1']\n",
    "    u2 = params['u2']\n",
    "    model = transit_model(x, epo, per, rprs, ars, b, u1, u2)\n",
    "    return np.sqrt((y - model)**2./yerr**2)\n",
    "\n",
    "def fit_transit_model(x, y, yerr, epo, per, rprs, dur, u1, u2, cap_b=True, fit_u=False):\n",
    "    '''\n",
    "    Fits transit model to data\n",
    "    INPUT:\n",
    "        x - array of times (days)\n",
    "        y - array of observed fluxes\n",
    "        epo - epoch (days)\n",
    "        per - orbital period (days)\n",
    "        rprs - planet radius divided by stellar radius\n",
    "        dur - transit duration (days)\n",
    "        u1 - first quadratic limb darkening parameter\n",
    "        u2 - second quadratic limb darkening parameter\n",
    "        cap_b - True if you want impact parameter within [0, 1];\n",
    "                False if you want impact parameter within [0, 1 + rprs]\n",
    "        fit_u - True if you want to fit limb darkening parameters; False if not\n",
    "    OUTPUT: lmfit best-fit results\n",
    "    '''    \n",
    "    best_chisqr = np.inf\n",
    "    # Due to degeneracies in fit parameters, try initial fits with b fixed\n",
    "    for b in [0.1, 0.3, 0.5, 0.7, 0.9]:\n",
    "        # Estimate ars from transit duration and impact parameter\n",
    "        nom = (1. + rprs)**2. - b**2. * (1. - np.sin(dur*np.pi/per)**2.)\n",
    "        denom = np.sin(dur*np.pi/per)**2.\n",
    "        ars = np.sqrt(nom/denom)\n",
    "        params = Parameters()\n",
    "        params.add('per', value=per, min=0)\n",
    "        params.add('epo', value=epo, min=0)\n",
    "        if cap_b:\n",
    "            params.add('b', value=b, min=0, max=1, vary=False)\n",
    "            params.add('rprs', value=rprs, min=0, max=1)\n",
    "        else:\n",
    "            params.add('b', value=b, min=0, vary=False)\n",
    "            params.add('delta', value=b-rprs, max=1)\n",
    "            params.add('rprs', expr='b - delta')\n",
    "        params.add('ars', value=ars, min=0)\n",
    "        if fit_u:\n",
    "            params.add('u1', value=u1)\n",
    "            params.add('u2', value=u2)\n",
    "        else:\n",
    "            params.add('u1', value=u1, vary=False)\n",
    "            params.add('u2', value=u2, vary=False)            \n",
    "        out = minimize(residual, params, args=(x, y, yerr))\n",
    "        chisqr = out.chisqr\n",
    "        if chisqr < best_chisqr:\n",
    "            best_b = b\n",
    "            best_out = out\n",
    "            best_chisqr = chisqr\n",
    "    # Now do a full fit, allowing b to vary\n",
    "    best_out.params['b'].vary = True\n",
    "    final_out = minimize(residual, best_out.params, args=(x, y, yerr))\n",
    "    return final_out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import lightkurve as lk\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "lcfs = lk.search_lightcurve(\"TIC 158588995\", mission=\"TESS\", author=\"QLP\").download_all()\n",
    "lc = lcfs.stitch()\n",
    "\n",
    "time = lc['time'].value\n",
    "qflag = lc['quality'].value\n",
    "flux = lc['kspsap_flux'].value\n",
    "\n",
    "good = (qflag == 0) & ~np.isnan(flux)\n",
    "\n",
    "time = time[good]\n",
    "flux = flux[good]\n",
    "\n",
    "mad = np.nanmedian(np.abs(flux - np.nanmedian(flux)))\n",
    "flux_err = 1.48*mad*np.ones(len(flux))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Normalized Flux')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEJCAYAAACDscAcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABFJ0lEQVR4nO2dfZwdVX3wv797dzeElxgIoRFiCJDQJk/BDayrqRKDlBXUygrtoxgFX8NWaLWKW/n0abWCoYhURBE3VShpRYpSecDypjyuoe4iBAORF4EQqxJIoUsppS0Ju/t7/piZy9nZmbkz987sfdnf9/OZz96dmTNzzsyZ8zu/l3OOqCqGYRiGkZZSozNgGIZhtBYmOAzDMIxMmOAwDMMwMmGCwzAMw8iECQ7DMAwjEyY4DMMwjEwUKjhE5EoReVpEHog5LiJymYhsF5FtInKMc+xMEXnM386MSHtj3HUNwzCM4iha4/hb4KSE4ycDy/1tPXAFgIgcAHwaeC3QC3xaRPYPEonIqcALxWTZMAzDSKKjyIur6mYRWZpwyinAJvVGId4lIvNF5JXAWuD7qvosgIh8H08AfUtE9gU+jidorkuTjwMPPFCXLk3KhmEYhhHm3nvv/TdVXRjeX6jgSMEhwK+d/5/w98XtBzgfuAT477Q3Wbp0KVu2bKkvp4ZhGLMMEfll1P6Wco6LSDdwhKp+N8W560Vki4hseeaZZ4rPnGEYxiyh0YJjJ/Aq5//F/r64/auBHhH5F+CfgSNFZDjqwqq6UVV7VLVn4cJpmpZhGIZRI40WHDcCZ/jRVa8D/kNVnwJuA/pEZH/fKd4H3KaqV6jqwaq6FHgD8Kiqrm1U5g3DMGYjhfo4RORbeI7uA0XkCbxIqU4AVf0acDPwFmA7ns/i/f6xZ0XkfOAe/1KfDRzlhmEYRmOR2TCtek9Pj5pz3DAMIxsicq+q9oT3N9pUZRiGYbQYJjiMwhkdHeXCCy9kdHS00VkxDCMHGj2Ow2hzRkdHOeGEE9izZw9dXV3ccccdrF69utHZMgyjDkzjMApleHiYPXv2MDExwZ49exgeHm50lgzDqBMTHEahrF27lq6uLsrlMl1dXaxdu7bRWTIMo07MVGUUyurVq7njjjsYHh5m7dq1ZqYyjDbABIdROKtXrzaBYRhthJmqDMMwjEyY4DAMwzAyYYLDMAzDyIQJDsMwDCMTJjgMwzCMTJjgKBCbasMwjHbEwnELwqbaMAyjXTGNoyBsqg3DMNoVExwFYVNtGO1ONVOsmWrbFzNVFYRNtWG0M9VMsWaqbW8K0zhE5EoReVpEHog5LiJymYhsF5FtInKMc+xMEXnM38709+0tIv8kIj8XkQdF5K+KynterF69mvPOO88+GKMumrHnXs0Ua6ba4miG+lCkxvG3wFeATTHHTwaW+9trgSuA14rIAXhrk/cACtwrIjcCu4EvqOoPRaQLuENETlbVWwosg2E0lLx77qOjo2za5H2SZ5xxRs3XCkyxQb7Cpthqx1uV0dHRKVaE8P9p09Vz/6bQ5FS1sA1YCjwQc2wION35/xHglcDpwFDcec7+LwEfTpOPY489VpuNkZER3bBhg46MjDQ6K0YdFP0eN2zYoOVyWQEtl8u6YcOGmq81MjKiXV1ditch0zlz5mTKd7is1crebnV8ZGRE586dq+VyWefOnatDQ0NT/k96DmnOS0Oe9SENwBaNaFMb6eM4BPi18/8T/r64/RVEZD7we3jCY8apt/fQNL2GGSKv3lazMRPvMc+e+/DwMC+99FLl/8CElCbPcWWNS9vodx7cf8GCBYyNjeWSj7D57frrr59mjou6R5TZrtZ241e/+hUdHV6z3UhNruWc4yLSAXwLuExVdySctx5YD7BkyZLc7p9HY1FrRZrJj7HtVOsCGB4eZvfu3UxOTrJ79+6aG4Qk8gqycBudQHgkNTzh95+lzo6OjrJ27VpeeuklOjs7C3kuSQR1Lng3pVKJOXPm1F33wkL8tNNO484776wq1OsR/q4A/NjHPsaePXsol8t8+MMfrsvUWC+NFBw7gVc5/y/29+0E1ob2Dzv/bwQeU9VLky6uqhv9c+np6dG6c+uTR++hloo0kw1wnvfKq7c1kyQJTffYggULmJycBGBycpIFCxYUItzrXc/EfZ/lcpn+/n4WLVoU2/BEvf8sdXbTpk3s2bMH8LSaTZs2zeg7D+qc+27yqHtRQvyoo46q+IyypEuD+x6CcnjWI68z3MjvqJGC40bgHBG5Fs85/h+q+pSI3AZsEJH9/fP6gPMAROQC4BXAhxqRYcjHdFBLRSqqAY5q6PK8V6OcpLU24OFG89JLL62YOoApx84880xKpVKlV7t169ZKr7CZtCv3fQL09vZy3nnnpTo/eP/nnXdey4SXB3Uu0DhEhI6OjlzqXpQQv/rqq9mzZw9XX3115DuvtS6G31tAVFnCQQ9B+sLeVZTjI48Nz5z0FPASnp/ig8AAMOAfF+By4HHgZ0CPk/YDwHZ/e7+/bzGeU+9h4D5/+1CavOTtHM/qJIxLl/WeYQdbvc7HOKddns684Hoz6SStJ/+u8xHQUqlUuc7AwEDlWKlU0t7eXp0zZ07k8ZlwXKYl6/Oo9/2PjIzonDlzVEQyO+DzYmRkRAcGBrSzs1NFRLu6ugrJRzVndT3PMkgrIpX6KCI6MDAw7Tw36KGzs3NKvayn3MQ4xwuNqmqWrcioqrQVI4/G2G2A87heUqVv5YiYLJEnUZ2Ajo6OykcYbOVyWQcGBnTu3LlaKpUqwqOrq0sHBgZyeydFETSkQV7TnF9vp6Ra+laKSItjaGhIOzs7tVQqRb7zevMQvLckQbBhw4YpwiUQMHmU2wRHQaRtfPOuxAMDA3VXjmZu6OohrXYWV/7e3t5pgsO9Tl9fX0V4tIrArfauZ7qhL7LuBflMGy6b9bphzbxUKmlHR4cODQ1NO79ao1/rvcPHTONoMcGR1twzNDSUm/oerijlcll7e3unVdy012rWhq4e818a7SxO+A4NDU0RGv39/akEzkySdQxFtQ5OUgPonlMulyta1tDQUM11J5yfgYGBXOphOJ/9/f2ptay01w3eebVn2tXVpSKiHR0dueShWv4CjXJoaCiTdpmECY4CifqIw43SwMBApSIl2VvTNJhhO7y7DQ4OFu5vyUIt98nD/FdN2wsL37AwHxoa0r6+vlhhXO14kSRpVEGjEe5xJj2rDRs2VDSooMeaZHJxtygTTbXe8cDAgPb391fy2NXVlVsPOZxPEal03Oqp7+51S6VS5d3HPdOBgYEpzynKLxGlAafNY1KHqVonIAsmOOogSkWt1tsLN0qu41REtL+/P7XpJOr6wXnuBx925lbLX5FOwyzlCZO2NxrX44vS9pJ6i1Hvo1rj1yiNI8pUFva9hH0z7nOJa2xcv06pVIp18oZt6dWefbh+h00qQc84z9Hx4XyWSiXt7OycUh+qaanh/90G2RWYcdfq7++fprWG85hWsMeVMao8aToBWTDBUSNJjVCcyh5ulAKV0f1oonpraf0gbq9tzZo1mZ1i1XpDeVGrkzpsbojrjUaZT0ZGRqZ8tK6GEWWfLpfL2tHRURHCaaaScLXJqCiXIgje+Zw5c6Y1Xv39/ZENuuubCT/jMK6T132W4TyE61u4Dic9m3C9E5FKfvL2RbjPKkmIRNUHtw64PfcooR1Vt6LKumbNmoqQDEfn9fX1TetYJtWpsPYTForVOgFZMMFRI2GTU19fX1WVPa5HEPWBu5XPDR+M84NE9Xw6Ojp0xYoVumbNmsTww+AeS5cuje0NRd0vzxDitOcFeV25cmWsMIwye4QjokqlUqy6HqQPh9/29vYm3jPcASiXy4WarOJ60YG5JJwft2EO8hVVJ8N28KQ6GBwL19/AtxaYycJ5Ca4R9dzC18/bdDo0NDRNC3O1dPfduj314NlG9dzTaLNBeTo7OyPfS+C8DkfnufVwzpw5U55rUqcnXJ6gExD4V+qpmyY4aiBc2cvlsg4ODsZ+qFFmgbCG4lYYoPJiXWEQbvjda4UrtNvrcPcHeXV78XH5jhNS1UIN0zw/12EXZx6I68VVa2iC51oul2N73EGvO+4jDKv2pVJpyv/hZxNn66/XJJBEOI+B3T5scovqzLgNo9tLddO4eY/qKQfPOuoe/f39kY1Y+JsI57O3t7dwE1/YXBTOr9tJCWtSg4ODUzohIqJ9fX06MjKiQ0ND2tvbW3G8x2nVUfd3y9/b2zvlmw+P13C/6eCdLV68WAcHByPbF/cbHRwcrGhb9WhxJjhqIOqD7erqiu1JRDXA4Q+2t7d3WiV1e9XhDy7cgLpmgKgKF9V4zJ07N7K3GCXwAsK9tawqb5IpKNz4RwnMcAMWbmziGsw4bTCsVQTmmMHBwUrvLOwzijIZuL25uEY6795zUpRXWICGn0XUGJOoZxSUM1w3g/vFPdew2SXc2CXZ8MPmySjBXisjIyORPp+o+0Z1qJYtWzatkxalGQCV+h1uoJMER3CtcEfSrXtR+13B5pY1rJHkZa4ywVED4Q+2WkPt2jEHBwcrQsKtHNUa+vAHF9WAArpixYpYARZ1vd7e3sjzo9TZcMWD7D3quIY9sGvHndPZ2an9/f2xH1Ogwsed09/fr0NDQ7pixYqK9pD0AQbXLZfL0wR4uMyuoCuXy9rd3T1NIOZtr1edai4NGoFw4xdlcgvb8l1TVNxzi3o+vb2905zvpVKpksYNM3c1u/7+fl2xYkVF0IXzHNbE3Wdfb8j6hg0bIssSNtts2DB98FwtW/Dtu413nIbv1rvu7m494IADph2LE/DBtmzZssSyJ9XjLJjgqIG4RjupIkTtL5fLU+zm1Ta3Zx3Xa4kzz8SZFAJfiKupBD2qsDobZb7JaieNssuHP96kc9xtzpw5qZ5zcG7QmKU5P24rlUra3d09ZXxMVKhneKxAloCAtM8xHKHnakmB0A/X1aVLl07Lx8jISOTgRvAER19fX+SxsGnEDfcNNPCweTUshFwfwYYNG6Y5g6PuW0/QQVzDvWLFiinnDQ0NJTbQWepL8I0NDg5qX19fZLmqdWLcupVUb12NI1zucJuxbt26mp8jJjiyE9Vo19o7WbFiReq0wQcTdrCFbaBRacO9l3nz5sVW1rCz1TW1BB9dPY7foIcbvr/byLoNUR49v0BI13ud8BZE1UQ5qV1TjRugkMc8TVFhw+GGrlQq6cqVK6ftC+pONd8ceA374OBg5LGw8AvnKSqQIEoIuWHD5XK5ohF1dXVF1tFaBUdQtwYHB/XQQw+NbXRdDTKPupdmy9KBDNfr5cuX67x586YJgsDnElg3oszeeWscLbcex0wxOjrKzTffPGVfqVSiVCoxPj4emaZcLk+byTLg4YcfTn3vXbt2VWa7dBfeCefFnWZZRFBVnn322SnnDQwM8OUvfzlyls0gbcCCBQsqv0WkMqvoUUcdlTrvLsGsnD/60Y+mlH/Xrl3TZqAN9v/TP/1TbJnT0NHRwcEHH1xz+jguvvhiPvnJT3LmmWeya9cubrnllko9CJ7r5OQkmzdvBl5+H/USnl3YvV/A5OQkDz300LR9wf0nJia45JJLYusmwEsvvcRdd90VWYfdehHOU7lcRkQYHx+nq6uLBQsWcOGFF9Ld3c3tt99eSdPZ2QlQmbEWvPq3fv16Vq1axTnnnFPZH5wfzPKahfBaHCIy7ZyLL76YI444grGxsSnTr2dl0aJFPP3006nTl0olDj744ClroqRlcnKSxx57DIDvfOc7nH322axevZqNGzdy1llnJaadmJjIfWp7ExwxDA8PT/uASqUSb3vb27jxxhsjK0u5XOb1r389d955Z12Nxg033MAtt9zCihUrpux3r3nkkUeyY8cOXnrppZfVxxD9/f1cdNFF9Pf3Mzw8zHXXXcd9991XOR58VKrK5OQkZ599NgBbt25lz549qCrj4+N1rVh2wgkn8OKLL047Fgiy3bt3c8455zAxMVERVvVwxBFHcPLJJ/O9732v0rCXy+VKGZMIhHEU27dv56yzzkJE6Orq4rLLLmNsbIxbb721IixcVJWJiYnc1oAIpszetWtX6rRunUgSGgFR5RARxsbGIvMUTNsNVBYb+qM/+qPKAk6Dg4MMDw9z8MEHMzg4CMDGjRsr1wme9djY2JT8rVy5kq9//es1PTd3cS2Y3jkK9n3kIx/hq1/9KuVyOfHcJNx3EdTdpDpWLpe56aabpqyJMm/ePO677z527NjB9u3bY9O6edu9ezdnnHEGp556Ktdcc820c9PkpV5McMSwdu3aab2vyclJFi1aREdHR2VxFZfx8XH22muvqr3Nzs5OVq1axT333BN73u7duys9jCgefvhhSqUSEF/h99lnH+DlNQTWrl3LG97whkqF6ujoYNWqVdx9992V/P/hH/7hlGvWs45BsJ5AVP6CHuvk5GSm3lepVOL000/nm9/8ZuTxhx56iIGBAY477jgAXnzxRXbv3s39998feb6I0NnZyQc+8AEAhoaGUFVEhN/6rd/iP//zP9m5c2elDKrK7t272bp1K2eccQY33HBDbF7zWgMCXl7zoR7Bunz58il1Ko1WJCL86le/YuPGjYyNjbFgwQK2bt3Krl27WLRo0ZRzb7nllikLOA0PD3PppZeyevXqigYdvt9Pf/pT5s2bN6WRS6r31XAX10piYmKCK664gvHx8bo1QxFhr7324tJLL+WWW26JrBOvfvWr2bZtW+VewZooo6OjPP/88/zgBz/IdM/t27fz+c9/ftr+sr864KpVq/jjP/7jiqZai/aWSJT9qt22vHwcgQ07zmcQ+AOiHNSLFi2aFmHi+i/ycNCFt1KpNG0Alzv4L4iMSbrG0qVLa54sbWRkZJpjG992vW7dutgyB07fuOgfVdV169bV/XxERFeuXFnx4YQduq5fI5x2+fLliY7OPEeUh30K9ZZ76dKlum7dutTXCoeAu5s7mC3qeOAsTgraCNcDN/KulmeV1gGd1+aGR8eFvbtlLJfL2t/fXxm4mpd/ZdGiRVMCNfIIC8ec49kJN6rVXnBnZ6cODQ1FxsK7H0p/f/80R2SpVNL58+fnXqnjBtMFW3d3d6rr1OLoHRkZiQy3TWr0y+VyxckX9/y6u7tzdWYGUWNuZFLgiO7r68vUyAZp85wDzH1/UQ1ttbzElTnvupY1D3H76wkqSBupl+cWdBDyitBKejZp0+Q1mwGNEBzAlcDTwAMxxwW4DG+lv23AMc6xM4HH/O1MZ/+xeCsGbvfTSrV81CI44gbeVXtxRWgOtVa2IJoiKVb9oIMOSn39rL3AKG1mYGBAly1bVrUcM/nhB1uejWmeI8kHBwenRXIFU8w0sr5FPTMRyRRBGLW9+tWvruvZJY1XyaO8CxYsmLIv0DiiFgCrdVuxYsW0SLnwc64WPZhHHaRBgmMNcAzxguMtwC14AuR1wE/8/QcAO/y/+/u/9/eP3e2fK37ak6vloxbBUYvKW1Rjt3Tp0tTaiDvQMDDDJM2bkyXPWXowIyMj064tIjoyMhIZ9tnoBjDvrR5zi0vcoDxgyhQfjRC0gO67777T6l8eDXY9IaTus0sz+DPLO426XjD79UybyKq1CfVOcKgaLzg872pBqOpm4NmEU04BNvl5vAuYLyKvBN4MfF9Vn1XVfwe+D5zkH5unqnf5hdoE9BeR97Vr11acz2nxspQ/S5Ys4RWveEXq87u7u1FVHnroIc466yxuuOGGilM17FzNkudwdE0SUY7Q4N4XXXQRfX19U46lifqJYuHChaxZs6amtEVSKpVycYxff/31scfGx8enRARlra9RiAjlcplyuZzq/BdeeGHK/0uWLKkrnDogCCGth7GxscRggiyBBqVSiVNOOSXS8b5nzx527dpVCTmeKZ577rnE45OTk9NCqfOiUMGRgkOAXzv/P+HvS9r/RMT+3Fm9ejWXX345nZ2ddYeI1svmzZv55S9/GXt85cqVld+Tk5NTQm4B/vEf/7HSMNcq3LI2hFFho6paaQzyija64IILppQ/II9GNA1RdaNUKvHVr341l7j50047LfW55XI5lRBdunRp5PMREU488UQ+/OEPZ8pjQGdnJ0uWLKl6TjAeqmiC8SZx91LVTN/2k08+GXm+qnLTTTfFju9qFFGh1HnRaMFRGCKyXkS2iMiWZ555pqZrrF+/nh/96EeceOKJOecuP1auXMlHP/rRxB7iqaeemroHGce5555bd0Ooqlx11VWMjo5W7S2l5fHHH+eMM86Y9kEfdNBBdZc5DWFBvGbNGs4991yuv/76KWMWamX9+vUMDg6mamjHx8c54IADqp47d+7cSr7d56aqdHd3A14ocTAeYPny5anyunz5cg444IDY46VSiQ9+8INccMEFrF+/PrHR7ujoqDuENBhv8ru/+7ux50xOTlIul1m8ePG0Y4ceemhlgOPk5CR33313bMdrYmKiZq25KFS1MI2j8IgmYCnxPo4h4HTn/0eAVwKnA0Ph8/xjP3f2Tzkvbqt3BcAkO3Ojt2CWzUWLFk075k5PUC3sNtjCNmvwIq+yEjfPVzAtRZ7TggwODkZGaq1bty52Ztei3kX4Oee1Tkcae32Qh2p2+6TjQdRW0jQ1SdPYJ0VxhWdGTnIm5+HjCIibKTecv6g8BGVK4wtJCtPPs56l9csU6eNotOB4K1Od43fry87xX+A5xvf3fx+g0c7xt1TLQz3rcQQTsjXK+ZjHVm2eojQVsJZQ3KhFlaJWKSvyYwrWAombwC/Pbe7cudMiYfr6+mqqe+FnuWLFitT5SGqowot4Za0H4bW203wX3d3dU9ZnD76rpFDwvIILgvslPZM0ArlaOUulkg4ODmp/f7/OmzdvyrHe3l4dGBjQhQsXFl4H3a2Vo6q+BTwFvITnj/ggMAAM+McFuBx4HC/EtsdJ+wG8kNvtwPud/T3AA36ar1BgOK67nkQRL3bx4sVVK2NUg59ViC1btqxuwZe0SmAc4WVEg4VwkhYfqresUVtXV1cuAwaTtr333jtycGjcLKZZ62HafATx+/39/dMm+AsGsNZaxkCjcqdoT9Ogukv/ZhE69Whr7uC3OO13/vz5kevjZA0HD86PSxfMLlzPs69lq2dW3ABsAGA20jZu9TRs1UYe52XOyaPRrMV0EL5v2unUFy5cWBkkmYfQg2zjVdJs+++/f6rz6u01ZxGyQZ3q7++PDYl1p9avpy4Ewqnaub29vVOmd49bejlqq2d2XHetj6Tw4KBzFsxo3NnZmaqBz/IMs45rycu6sWzZstbUOJplq1fjqMdGWat6GhejnWUK6EB9jhoAmLVyZrWXhv1CYdND0ijbYOR1eAGsPD4mtzwrVqxINN0EU6THHat2jzxMBeFVEsPboYcemosg6O7unmZiqVa2amavYCaFuDW6q5lPaxEcwZojcdOkRH1XwUqJwfr1eZiuqm1x7UIwFUle96t3BgNMcGTHXTMibgqMah9knrb17u7uSmNaKpV06dKllVXYwpVtzZo1lXlrwj6FWiph1kYwSltyF0WqZuZw1X4R0UMOOaSqaS/PLbDnxwm3NNOe5DXtSNwCTCKig4ODqR3nUZu7+mPWIJBw+cOCYOnSpZX8Bw1zeJnYOMFci4YbNbVOOI9hLVj8RcwCgdaoGQvcbyONMz/LVs+cadh6HNkJZpUFOOqoo9i0aRPf+MY3Ug9wmpiYmLImQT0sX76cgw46iG3btjE5OUmpVOLII4/k5JNPZmxsjEWLFk2Z7dSd2r2zs3PKwKVgfxayTtEctSbG3Xffzd13383jjz/OlVdeOS0fbv5VtTI7saqyc+fOzHmuFRFhzpw5nHbaadx55528+OKL0/IajJVJmmE2j2nVw/dzUVU+//nPs27dupreqYhw+eWXs379ejZu3MjFF1+cKb17z6VLl6KqU8YbzZ07t/I7mN336quvrqy/snXr1kq9EhGOO+44fvzjHzM5OUlHR/amKZiNOS6PpVKJZ555Zto7O/XUU7n++uvZvXt3Tc8RYP78+TWHmPf09PCTn/xkyr7f/M3fzLSGz4wTJU3abas3HNdlZGSk0AidpInfqvlEispTcP0spip3mpNw3pYtWxa7HnleZarHvOjOdlrP+85rzfG49bPrfU6lUqlinsnj2uFzo5bcdddCD8yRgRnLXU62lmV3k6KnAvNntfDs8BLLUdeRhNmba9lcTSNYDTMvjaOWiEgXzFSVnbhpiaup9LV+yFGOuWBa83oqUuD4q1bBxF/uNFgzOdxg1rLueLCEp3udYJrtqEkZ036MgY8iy3Mol8uVhsqd5ym8DQwMTHnvSZFNcT4wV/jUS9ySrmme0X777Rd7PGi4w/sPOeSQiolzw4YNsf6cqDoVNrsEzy+YPtwdI+L6FoL7uf6QWp5f2HcW3HPFihUVgRUnGILIv8A0HTU31YoVKyomt6SxSkEZws/OXZs8vMxC1DexePHiqt9+XEcs+NbqARMc2XCd4+FKnDTbrPjTcUc1gIFfIpy2o6MjtjIGH1QtDQegK1euTDUOoLOzc8pc/qrTB/HVYysNxlLE9a4CW3vcxxg822Dr6urSoaEh7erqSpyVOPj4A6EYNIYjIyOVtZqXL18+7X247z3Ia7i3GrzrKP9DXmtxqE5fvzsIegimfE/yE0TlLdA03B5+XGMTN1Av8K9FXTusKQS+h+A9BcIj7CwP/q93DYngfcX59YL1MMLfYeDHc/1vQUMffq7Bdxnl4O/u7q4IT1douvXPrbuBJhTVpgTrnSR1/II15cN1s16hoWqCI/MDc9XrsNpcLZx0aGio6tz8QRhgeJEk9yMLnKtJgiqpkXXXNQgLn3DeXBNC0GDmITiSFpMJO01Vp/eu3ciwKFNG0iDNYBR0sN9NEwitpHE6QQPrNmzr1q2rNCbhyC93q2dNiTDhd9fX11dpfIKIq2Adk2BxILchDi8YFjxPt/EKyhQ4i4O8x9W9oBF0G87w6PCAqJmmg45C0ndWD3EdkKB8g4ODU/IUNLRhB3vceBPXqR4XNh9o7729vdM00HD+3Ly4kZOuVpbUnoS/gbyeZZzgMOd4DMEEacHSi2vXrmV0dLSyznKwDvRVV101ZXnUUqlUmVgszqEsIhx++OF88pOfZP369ZX9wfW//OUvMzY2xtq1ayuO1XK5nGkStRNPPJHPfOYzlfTBfa6//npOO+20irP/qquuYnx8HBFhYmKCycnJyrKfZ5xxRqV8tSw/Gaw5HqS/4447KvkJH1u1ahXDw8N873vfm3KNnp4eLrrookqaq6++mt27dyMiPPfcc1PWvb766qunOLKDuZCC/8GbUTRYB919Z1HrjU9OTrJr164p66Nfd911lTRvectbuOmmmyrBCgcffHBlmdl61moPs379eh5//HG+8IUvMDk5ye23384PfvADSqVSJXgAvMkw58yZU1kPPag/wbsGWLVqFR/72MfYs2dPZS324FkFjULw/oPlhjs7O6c5ncfHxxkbG2N4eLiy3ni4zgZELeeqqpXzw99ZUXR0dPChD32IVatWcfbZZ08LGPnSl77E888/P225aLf+uOfv2bOHsbExLr30Uo4//nh279495Zzdu3dPWd715JNPjq0PbpDAG97whilBAqtWrWJsbIxPfOITfPGLX5wWnFMul6d8AzPxLBuuDczEloePIzyoyA0tjDK5hHsuQY8+6FlEzd2TZN8Neo7ugCVieh9ZQhnd0Mio+9ez/GRSbzJsDnCfjbuF/SrBc3DTuialvr6+2LEfgVMzqvcbaBLh+wcaR9hkEaWNDA4OVsqRl2M86lm65Yl6/2vWrJmmyUVdx00fpXEGuKY6V5tJKl945LZ7r7BWU08dS7p/eCG2QEvq6+uL1SLcd5rkJHffcZx1IZzWnX4mLk1nZ+c0DTGsQbqLVIVX+gu3WbZ0bAMEh0v4ww1X/nDDHthAw/HrbsPmNqZpVHa3IsRVvHBFykLeH3CcMIwyB0Q12lHTnESZpNznFRYsUR98VAMcZb4AKuMkogR10EEIm43csRF5ER4IWC36J9jCJrM4M6sbFJGUd7e8SSbIuJHbnZ2dUwRalLkyz2fmfm9R7z7qWYXfqTueKFhD3W2co0yd5XI588wJwWwRcdGFYVNrlJk7qqNbTyfGBEeNxFWipIY/SrC41wu/0JGRkYrjNW1vNWw3DnrAcemK6NWlIeq+4V5vf3//lGcb9wzCAid8rtu4hkfYu5EuUY7IwJYctlcHDV2UsAnm3ooqU57Ocbf8Qc/fjUKqFqYd7oREDRg89NBDU9e/ao1SnFYTfi7h95mnXyic16gZF4JIK7fDFW6U3e8/qrGOWyk0aBsCDccVNnHTyAQj7cOCNm6urySLRb2hzQE1Cw5gr4h9B1ZL10xbPaaqcMSH61SNenFxgsXFjTCK6n27PdkkQeD2dJKmBMmr95EXUflxNYUg6iXcmwqbofr7+6c8o/B4gcCxGY4Yi+q9xzm6XZNU2AQWFloz1QhGRXzF9aSj8pEUuRauS1H1sJoTNvx84/ISdrxnHStU7Vm52nmUxrhy5cppTui4d5ikPbsBClHmvqh2JFyf4kzca9asmaLdxVkmwvvDZtQZ1zjwZq19nfP/acCj1dI101ar4Ih7SXGNevARJ9mBo3oGUb2yNC/dNY8lnVdU5Epaop5XeF+SxhbV0LvH4vxQbmMRZ1cPPkq38XWjjML3CYR+lLmxiKgWlzSNRiBQ48w/Ubb/qJ5vnLkjrYCMEvRhLawoYZvU+3bLGPYlhDVRV5AltQVRGkqchh2OBgyb/KI0mHAYbxqNIy4vWalHcBwF3ANcDHwTuBVYXC1dM215aRxpX0DSC4v6yMMfjzvgr9rkgmkqRyM1jrT3TtLYwppEYCJK+liiBk1m0cqSbPhxArvo51yt15vmvlH1L2yjD+L/oxq8LCa5NPkKN7x5kNT7DvyPYWEShL+7wqSaxpG24+Z2fKr5v8KWBFfQVaubeQiKMDULDi8t/cB/Ak8Cy9KkaaYtDx9HXi8jqpcVtp+G4/bTOFrTmLaq9fqLIIu2E6exxTVA1cbahD/AahM1ZhHCcY1A0c80SdtNc9+oZ+marsJmqijBnEU4ZslXXs8tTe87rqPijmRPesbh+lWtg5dWyATnukLNDR6ZaYtBPRrHN4Bh4DDgzcDPgbOrpWumLc+5qtJQ7SOoZtJw1dVqFTK4X9aeblKaoj/iNGnSCLlq13b9JnlFOjXa7JcH4QYwybw6Ex2OIjS1tJ2AaqblOMImpWqdkjituVre8pqKpVbqERwfw1llD3gF8I1q6Zppm0nBkVY9Tzon64dUS2OWZLNtxEdc1LVboZFrFG5ZomYxmCkaLYxrqSNZzE/h88N+uqLymAd1mapq3YCTgEfwln/9VMTxQ4E7gG2+VrPYOXYR3hKxDwDvdPafAPwUuA/45zSms5kUHGk/gjwbvDw1jkZ/xK1Aoz7iPHDz3izvulWFcdZ6MDIydZGpVvi+6tE4fgHsCG8p0pXx1gU/HOgC7gdWhs75NnCm//tNwN/5v98KfB/oAPbBc87P8489Cqzwf38E+NtqeWk2jaOo+9bSa8pq/jGm0kpCJPxuG2kCicpbqzzHWhkZKT5kO2/iBEeauap6nN97AX8AHJAiXS+wXVV3AIjItcApwEPOOSuBj/u/fwjc4OzfrKrjwLiIbMPTXq7zH/o8/7xX4Dnsm4bVq1dzxx13VOZQymOuoqJwF6py97VK/htN0lxczUiw0NHExERlnqVmeddRdbHdGB4ersyDJSK8//3vb90yR0mTahtwb4pzfh/4uvP/e4GvhM65Bvio//tUPKGwAOgDfgzsDRyIp+V8wj/vOGAMeAJPCM2rlpeZdo7PNO2qJTR7L7RZTD1padd60irEPf9mrufUqnGIyDHOvyU8DSSvWXXPBb4iIu8DNgM7gQlVvV1EXgOMAM8Ao0AwZeWfAG9R1Z+IyCeBvwY+FJHv9cB6gCVLluSU3eYk3JPMa1bWRtIKvfmZnNk1D0ybbCxRz78V6nkUaQTAJc7vceBfgP+dIt1O4FXO/4v9fRVU9Uk8TQMR2Rc4TVWf8499Dvicf+wa4FERWQi8WlWDBXr/AW9A4jRUdSOwEaCnp0dT5LdlabUGLA2tIAxbsSGeDSahvHGXU6j32YWffyvU8yiqCg5VPb7Ga98DLBeRw/AExruAd7sniMiBwLOqOgmcB1zp7y8D81V1TESOBo4GbveTvUJEjlTVR4ETgSZe0T0/kipvKzZg1WgVYWgNcXtTtEbQKvU8TKzgEJGPxx0DUNW/rnJ8XETOAW7Di7C6UlUfFJHP4tnNbgTWAhf64WmbgbP95J3Anf5CPM8D71HPUY6IfBi4XkQmgX8HPlC1lC1Omsrbig3YbBOGRutRtEbQqvU8SePYr96Lq+rNwM2hfX/h/P4O8J2IdC/iRVZFXfO7wHfrzVsr0arqbBLtKgyN9mImNIJWrOexgkNV/3ImM2LE06rqbBLtKAyN9qNVNYKiSTJV3a6qff7v81T1wpnLluHSjpW3HYWh0Z60okZQNOKF6kYcENmqqqv83z9V1WMiT2wBenp6dMuWLY3OhhEiz2gVwzDyR0TuVdWe8P4kH0dbh7C2O63QKFtPzjBakyTBcbiI3AiI87uCqr690JwZNdOqg4oMw2gNkgTHKc7vLxSdESM/zPFsGEaRJEVV/WgmM2LkhzmeDcMokrzmnDKaiFqisFrBJ2IYrUK7f08mONqULI5n84kYRn7Mhu+p1OgMGI0nyidizD5GR0e58MILGR0dbXRWWprZ8D0lDQC8iYSQXIuqah/MJ2LMhl7yTDEbvqckU1UQSXUqsAj4e///04F/LTJTxszSjiPTjWxYJF5+zIbvKXbkeOUEkS3hkYNR+5oZGzluGMmYxmFEUcvI8YB9RORwfXnt8MOAffLOoGEY+ZIlsmc29JLbhWaI2EojOP4EGBaRHXijyA8Fzio0V4Zh1EUtGoRNAdP8NItmWDWqSlVvBZYDHwX+GPhNVb2t6IwZhlE7syGyZzbSLO+1quAQkb2BTwLnqOr9wBIReVvhOTMMYxppQ2aDyJ5yudy2kT2zkWZ5r2lMVVcB9wKBPrQT+DbwvWoJReQk4Et4S8d+XVX/KnT8ULx1xhcCz+ItEfuEf+wi4K3+qeer6j/4+wW4APgDYAK4QlUvS1GOlqYZ7JpGY8lipjCfRXuyevVqLr30Uq6//npOO+20xr1XVU3c8NYHB9jq7Ls/Rboy8DhwONAF3A+sDJ3zbeBM//ebgL/zf78V+D6eYNsHuAeY5x97P7AJKPn/H1QtL8cee6y2MiMjIzp37lwtl8s6d+5cHRkZaXSWjAawYcMGLZfLCmi5XNYNGzY0OkvGDDPTbUHQ/oe3NCPH94jIXPzBgCJyBLA7RbpeYLuq7lDVPcC1TJ1xF7x1xf+f//uHzvGVwGZVHVfV/wK2ASf5x/4Q+KyqTgKo6tMp8tLSzIRd00YNNz/NYqYwGkez+DjSmKo+A9wKvEpEvgm8HnhfinSHAL92/n8CeG3onPvxBhh+CXgHsJ+ILPD3f1pELgH2Bo4HHvLTHAG8U0TeATwD/LGqPpYiPy1L0SNRmyVSw0jGzE9Gs4xKryo4VPV2EbkXeB1eOO5HVfXfcrr/ucBXROR9wGY8/8mEf8/XACN4wmEUz58BMAd4UVV7RORUPB/JceELi8h6YD3AkiVLcspuYyi6wbBRw62DhczObpql85Bm5PgdwCWqerOzb6Oqrq+SbjXwGVV9s///eQCqemHM+fsCP1fVxRHHrgH+XlVvFpGfAyer6i98R/lzqvqKpLzYyPFkTOMwDCOKekaOHwb8qYi8RlX/0t+XZrqRe4Dl/kjzncC7gHeHMnUg8KzvrzgPT3tARMrAfFUdE5GjgaOB2/1kN+CZrn4BvBF4NEVejASapRdjGEZrkEZwPAecAFzmz5j7njQXVtVxETkHuA0vwupKVX1QRD6L56m/EVgLXCgiimeqOttP3gnc6SkUPI8XpjvuH/sr4Jsi8ifAC8CH0uTHSMZMIIZhpCWNqWqrqq7yf78P+ASwf5RJqVkxU5VhGM1Ms47TqsdU9bXgh6r+rYj8jJc1A8MwDKMOWtHHGDuOQ0Tm+T+/LSIHBBueb+HcGcmdYRhGm9MsYzOykKRxXAO8DW+6EcULxQ1QvBHhhmEYRh00y9iMLMQKDlV9m//3sJnLjmEYxuyiFaMak9YcPyYpoar+NP/sGIZhzD5aLaoxyVR1ScIxxZuU0DAMw5hlJJmqjp/JjBiGkT/NGuZptDZpwnERkd/Gm7F2r2Cfqm4qKlOGYdRPK4Z5Gq1BmhUAPw182d+OBz4PvL3gfBmGEUGW6e+zhnna1PpGWtJoHL8PvBpvIaf3i8hvAH9fbLYMwwiTVYPIEuZp2omRhTQLOf2PPwnhuD8o8GngVcVmyzCMMFk1iCDM8/zzz68qCFpxEJrRONJoHFtEZD7wN3iDAV/AWx/DMIwZpJaBYmnDPFtxEJrROKpOcjjlZJGleGt/byssRwVgkxw2FovsyY8in6W9JyNM3CSHqQSHvybGUhwNRVX/Mc8MFokJjsZhtnPDaF1qnh1XRK7EW0jpQWDS361AywgOo3HUsiyt9XwNo7lJ4+N4naquLDwnRluS1XbejBqKCTLDmEqaqKpREalJcIjISSLyiIhsF5FPRRw/VETuEJFtIjIsIoudYxeJyAP+9s6ItJeJyAu15MuYObJE9kDzRfcEguzP//zPOeGEE6qOcbCxEMZsII3GsQlPeOwCduNNr66qenRSIn/d8MuBE4EngHtE5EZVfcg57QvAJlW9WkTeBFwIvFdE3gocA3QDc4BhEblFVZ/3r90D7J+hnEYDyTKBW7NF92QxtTWjtmQYRZBGcHwDeC/wM172caShF9iuqjsARORa4BTAFRwrgY/7v38I3ODs3+yvMz4uItuAk4DrfIF0MfBu4B0Z8mO0AM02xXQWQVaLP8cwWpE0guMZVb2xhmsfAvza+f8J4LWhc+4HTgW+hCcE9hORBf7+T4vIJcDeeFOdBALnHOBGVX1KRDDaj2aaYjqLIGs2bckwiiKN4NgqItcAN+GZqoDcwnHPBb4iIu8DNgM7gQlVvV1EXgOMAM/gDTicEJGDgT8A1la7sIisB9YDLFmyJIesGrOVtIKs2bQlwyiKquM4ROSqiN2qqh+okm418BlVfbP//3l+wgtjzt8X+LmqLo44dg3e/FiCZzp70T+0BNihqsuS8mLjOAzDMLJT0zgO358wpqrn1nDPe4DlInIYnibxLjy/hHv9A4Fn/bmwzgOudO47X1XH/MGHRwO3+z6PRU76F6oJDcMwDCNfEgWHqk6IyOtrubCqjovIOcBtQBm4UlUfFJHPAlt8v8la4EIRUTxT1dl+8k7gTt+H8TzwHl9oGIZhGA0mjanqCjxH97eB/wr225QjhjG7sYGR7U/NU47grfo3xtQ1xm3KEcOYxdQyZsUETftQVXCo6vtnIiOGYbQOWces2ODI9iLN0rGLReS7IvK0v13vTg1iGMbsIxizUi6XU41ZabapZIz6SGOqugq4Bm/8BMB7/H0nFpUpwygSM5nUT9YxKzY4sr1I4xy/T1W7q+1rZsw5bgSYyaRxmMBuPepxjo+JyHuAb/n/n47nLDeMlsPmk2oczTSVTNG0u5BMIzg+AHwZ+CJeNNUIYA5zoyUxk4lRNLNBq00TVfVL4O0zkBfDKBybT8oomtmg1cYKDhH5i4R0qqrnF5Afw8hMVrPAbDKZGDPPbNBqkzSO/4rYtw/wQWABYIIjRLvbNZuR2WAWMFqL2aDVxgoOVb0k+C0i+wEfxfNtXAtcEpdutmINWGOYDWYBo/Vod602cQCgiBwgIhcA2/CEzDGq+qeq+vSM5K6FqGWAk61PXT9ZB6IZhlE/ST6Oi/FW59sIHKWqL8xYrlqQrHZN01DyYTaYBQyj2UjycXwCb8W//wP8mbNMq+A5x+cVnLeWImsDZiaW/Gh3s4BhNBtJPo6q81gZU8nSgM2GyAvDMNqTNAMAjQIwE4thGK2KCY4GYiYWwzBakULNUSJykog8IiLbReRTEccPFZE7RGSbiAy707WLyEUi8oC/vdPZ/03/mg+IyJUi0llkGQzDMIypFCY4RKQMXA6cDKwETheRlaHTvgBsUtWjgc8CF/pp3wocA3QDrwXOFZHAGf9N4LeAo4C5wIeKKoNhGAZY6HyYIk1VvcB2Vd0BICLXAqcADznnrAQ+7v/+IXCDs3+zqo4D4yKyDTgJuE5Vbw4Si8jdgC0qZRhGYVjo/HSKNFUdAvza+f8Jf5/L/XhjRQDeAewnIgv8/SeJyN4iciBwPPAqN6FvonovcGsBeTcMwwBs9cIoGh1yey7wRhHZCrwR2AlMqOrtwM14U7h/CxgFJkJpv4qnldwZdWERWS8iW0RkyzPPPFNYAQzDaG9sdoLpFGmq2slULWGxv6+Cqj6Jr3GIyL7Aaar6nH/sc8Dn/GPXAI8G6UTk08BC4Ky4m6vqRrxR7/T09CQvc2jMKmwySiMLFjo/nSIFxz3AchE5DE9gvAt4t3uCb4Z6VlUngfOAK/39ZWC+qo6JyNHA0cDt/rEPAW8GTvDTGUZqzF5t1IKFzk+lMFOV79g+B7gNeBjPsf2giHxWRIKFodYCj4jIo8Bv4GsYQCdwp4g8hKc1vMe/HsDX/HNHReS+KuuGGMYUzF5tGPVT6ABAPwLq5tC+v3B+fwf4TkS6F/Eiq6KuaYMWjZqxqV4Mo36sETZmFWavNoz6McFhzDrMXm0Y9dHocFzDMAyjxTDBYRgJ2FQThjEdM1UZRgwWutu+2Fie+jDBYRgx2CqN7Yl1COrHTFVG4bSqucemmmhPbCxP/ZjGYRRKK/fuLHS3PbGxPPVjgsMolJkw9xRpr7bQ3fbDOgT1Y4LDKJSie3etrNEYjcM6BPVhgsMolKJ7d+bANoyZxwRHjliIXzRF9u7MXm3Ugn2r9WGCIyeazWQyWz4Ms1cbWWm2bzXIUyvVYRMcOdFMJpNm/DCKxOzVRhaa6VuF1vxebRxHTmSN+S9ybIPFqRtGPM02PqcVv1fTOHIii8mk6B5GLXb/VlOVDaNWms282ZJ+OlUtbANOAh4BtgOfijh+KHAHsA0YBhY7xy4CHvC3dzr7DwN+4l/zH4Cuavk49thjtZnYsGGDlstlBbRcLuuGDRtyv8fIyIhu2LBBR0ZGUp07d+5cLZfLOnfu3FRpDMPIjyzf60wCbNGINrUwjcNfN/xy4ETgCeAeEblRVR9yTvsCsElVrxaRNwEXAu8VkbcCxwDdwBxgWERuUdXnfYHyRVW9VkS+BnwQuKKochTBTPQwstj9m83m28qY5mbUQqv56Yo0VfUC21V1B4CIXAucAriCYyXwcf/3D4EbnP2b1VtnfFxEtgEnici3gTcB7/bPuxr4DC0mOExVbk9a0clpGLVQpOA4BPi18/8TwGtD59wPnAp8CXgHsJ+ILPD3f1pELgH2Bo7HEzgLgOd8gRJc85DCSlAgzdTDaDZB1qqY5mbMFhrtHD8X+IqIvA/YDOwEJlT1dhF5DTACPAOMAhNZLiwi64H1AEuWLMkzz21JMwmyVsU0N2O2UKTg2Am8yvl/sb+vgqo+iadxICL7Aqep6nP+sc8Bn/OPXQM8CowB80Wkw9c6pl3TufZGYCNAT0+P5lYqw4jBNDdjtlCk4LgHWC4ih+E17u/iZd8EACJyIPCsqk4C5wFX+vvLwHxVHRORo4GjgdtVVUXkh8DvA9cCZwL/t8AyGEYmTHMzZgOFDQD0NYJzgNuAh4HrVPVBEfmsiLzdP20t8IiIPAr8Br6GAXQCd4rIQ3haw3scv8afAh8Xke14Po9vFFUGozG06sJPRvtidXIq4oXqtjc9PT26ZcuWRmfDSIFFJhnNxmyukyJyr6r2hPfblCNGU9GK0y8Y7Y3VyemY4DCaimabR8gwrE5Op9HhuIYxBYtMyhcbyV4/VienYz4Ow2hTZrNt3sgH83EYxizDbPNGUZjgMIwWIktYqNnmjaIwH4dhtAhZTU9mmzeKwgSHYbQItUyiaCPZjSIwU5VhtAitbnqy0dftg2kchtEitLLpySK82gsTHIbRQrSq6cnWKmkvzFRlGEbhtLqZzZiKaRyGYRROK5vZjOmY4DAMY0ZoVTObMR0zVRmGYRiZMMFhGIZhZMIEh2EYhpGJQgWHiJwkIo+IyHYR+VTE8UNF5A4R2SYiwyKy2Dn2eRF5UEQeFpHLRET8/aeLyM/8NLf665YbhmEYM0RhgkNEysDlwMnASuB0EVkZOu0LwCZVPRr4LHChn/Z3gNcDRwO/DbwGeKOIdABfAo7302zDW9fcMIw6sZHdRlqKjKrqBbar6g4AEbkWOAV4yDlnJfBx//cPgRv83wrsBXQBAnQC/+r/FmAfERkD5gHbCyyDYcwKbGS3kYUiTVWHAL92/n/C3+dyP3Cq//sdwH4iskBVR/EEyVP+dpuqPqyqLwF/CPwMeBJP8HyjuCIYxuzA1u5oHZpBM2z0OI5zga+IyPuAzcBOYEJElgErgMDn8X0ROQ64C09wrAJ2AF8GzgMuCF9YRNYD6wGWLFlSbCkMo8UJRnYHGoeN7G5OmkUzLFLj2Am8yvl/sb+vgqo+qaqnquoq4M/8fc/haR93qeoLqvoCcAuwGuj2z3lcvTVvrwN+J+rmqrpRVXtUtWfhwoW5Fsww2o1gZPf5559vZqomplk0wyI1jnuA5SJyGJ7AeBfwbvcEPyLqWVWdxNMcrvQP/Qr4sIhciOfTeCNwqX+dlSKyUFWfAU4EHi6wDIYxa7CR3c1Ps2iGhQkOVR0XkXOA24AycKWqPiginwW2qOqNwFrgQhFRPFPV2X7y7wBvwvNlKHCrqt4EICJ/CWwWkZeAXwLvK6oMhmEYzUSzzPklnsWnvenp6dEtW7Y0OhuGYRgthYjcq6o94f02ctwwDMPIhAkOwzAMIxMmOAzDMIxMmOAwDMMwMmGCwzAMw8iECQ7DMAwjE7MiHFdEnsEb89GKHAj8W6MzMYNYeduX2VRWaI/yHqqq06bemBWCo5URkS1RcdTtipW3fZlNZYX2Lq+ZqgzDMIxMmOAwDMMwMmGCo/nZ2OgMzDBW3vZlNpUV2ri85uMwDMMwMmEah2EYhpEJExwNQESuFJGnReQBZ99nRGSniNznb2/x9y8Vkf9x9n/NSXOsiPxMRLaLyGUiIo0oTxJZyuofO1pERkXkQb9se/n7m76skPndrnP23ScikyLS7R9rx/J2isjVfrkeFpHznDQnicgjfnk/1YiyVCNjWbtE5Cq/rPeLyFonTUu820RU1bYZ3oA1wDHAA86+zwDnRpy71D0vdOxu4HV4i13dApzc6LLVWdYOYBvwav//BUC5VcqatbyhdEcBj7fSu63h/b4buNb/vTfwL379LgOPA4cDXcD9wMpGl63Osp4NXOX/Pgi4Fyi10rtN2kzjaACquhl4tp5riMgrgXmqepd6tXET0J9D9nIlY1n7gG2qer+fdkxVJ1qlrFDXuz0duBZa591C5vIqsI+IdABzgT3A80AvsF1Vd6jqHrzncEoR+a2HjGVdCfw/P93TwHNATyu92yRMcDQX54jINl8l3t/Zf5iIbBWRH4nIcf6+Q4AnnHOe8Pe1ClFlPRJQEblNRH4qIoP+/lYvK8S/24B3At/yf7dreb8D/BfwFN7y0F9Q1WfxyvZrJ22rlTeqrPcDbxeRDvGWzz4WeBXt8W5NcDQRVwBHAN14H9Yl/v6ngCWqugr4OHCNiMxrSA7zI66sHcAbgHX+33eIyAmNyGDOxJUXABF5LfDfqvrA9KQtSVx5e4EJ4GDgMOATInJ4IzKYI3FlvRJPKGwBLgVG8MreFpjgaBJU9V9VdUJVJ4G/wfvIUNXdqjrm/74XzxZ8JLATWOxcYrG/r+mJKyveh7ZZVf9NVf8buBnPptyyZYXE8ga8i5e1DWjf8r4buFVVX/LNNz8GevDK9irnEi1T3oTvdlxV/0RVu1X1FGA+8Cgt/m4DTHA0Cb7tM+AdwAP+/oUiUvZ/Hw4sB3ao6lPA8yLyOj8q4wzg/85wtmsirqzAbcBRIrK3bwd/I/BQK5cVEsuLiJSA/43v3wBo4/L+CniTf84+eA7inwP3AMtF5DAR6cITpDfOXI5rJ+G73dsvIyJyIjCuqi1flys02js/Gze83uVTwEt4vewPAn8H/AwvquhG4JX+uacBDwL3AT8Ffs+5Tg9eRX0c+Ar+gM5m2rKU1T//PX55HwA+30plrbG8a4G7Iq7TduUF9gW+7b/fh4BPOtd5C16P/HHgzxpdrhzKuhR4BHgY+AHeLLMt9W6TNhs5bhiGYWTCTFWGYRhGJkxwGIZhGJkwwWEYhmFkwgSHYRiGkQkTHIZhGEYmTHAYLYOILHBmId3lzEr6goh8taB7fkxEzvB//62I/MK/589F5NPOecMi0uP//hd/9tOfichDInKB+LP8ZrjviSJyr3+Ne0XkTRHn3JhyptbMs/CKyAEi8n0Recz/u7+/X/zztvvTbBzj718oIrdmfLxGi2KCw2gZ1Jv0sFtVu4GvAV/0/99XVT+S9/38QYgfAK5xdn/Sv383cKY/D1EUx6vqUXgjiQ8HhjLe/t/wxuwcBZyJN17AzdupwAsR6YJn0q2qNwOo6jed5/Ze4Beqep9//hXAh/EGli4HTvL3fwq4Q1WXA3f4/wOc7Jy73k+Pqj4DPCUir89YTqMFMcFhtDwislZEvuf//ox4az7cKSK/FJFTReTzfq/6VhHp9M87VrxJI+8Vb1LFV0Zc+k3AT1V1POJYoEH8V1LeVPUFYADoF5ED0pZJVbeq6pP+vw8Cc0Vkjp/3ffHmLbsg7fUc0s7Cewpwtf/76tD+TepxFzDfeXY34M0zZrQ5JjiMduQIvEb/7cDfAz/0e+7/A7zVFx5fBn5fVY/Fm5DucxHXeT3eOgouF4vIfXgjh69Vb86lRFT1eeAXeL30WjgNT4Dt9v8/H28yvf+OODevWXh/Q73pMQB2Ab/hpImbyXYLcBxG22OCw2hHblHVl/CmgigDge39Z3hTQfwm8NvA930h8H+YOvFcwCuBZ0L7AlPVIuAEEfmdlHmqaZU3EflfwEXAWf7/3cARqvrdiNMLmYXX10bSTDHxNN7Mt0abY4LDaEd2A6g3Y+lL+vK8OpN4U7cL8KDjCzhKVfsirvM/vGySmoJvghrGm/49ERHZD09gPRraf7bjsJ7W4IrIYuC7wBmq+ri/ezXegkD/AvwzcKSIDPt5ynMW3n8NTFD+36edNHEz2e6F98yMNscEhzEbeQRYKCKrobIW9v+KOO9hYFnUBXzH+WvxJqqLxfdHfBW4QVX/3T2mqpc7wuvJULr5wD8Bn1LVHztprlDVg1V1KZ7QelRV1/pp8pyF90Y8pzz+X3f/GX501euA/3BMWke69zTaFxMcxqxDveVJfx+4SETux5t5OMrkdAveOtMugY9jG57p6x/9/R34mo7PD/1Q2bvxphM/K2M2z8ETWn/haCUHVUkTBAFsA44H/sQ5tgb4taruCKX5CPB1YDueELzF3/9XwIki8hjwu/7/4K2RssM//2/89AHH4wk7o82x2XENIwER+S4wqKqPJZwzB68h/W1V/Y8Zy1yTISKbgVPCmpXRfpjGYRjJfArPSR6JP+jvPuCrs1xoLAT+2oTG7MA0DsMwDCMTpnEYhmEYmTDBYRiGYWTCBIdhGIaRCRMchmEYRiZMcBiGYRiZMMFhGIZhZOL/A6bSkeoUKQ07AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(time, flux, \"k.\")\n",
    "plt.xlabel(\"Time (BJD - 2457000)\")\n",
    "plt.ylabel(\"Normalized Flux\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get initial model guesses from BLS\n",
    "guess_epo = 1641.404552\n",
    "guess_per = 1.977143\n",
    "guess_dep = 0.0144563\n",
    "guess_dur_hrs = 1.0553\n",
    "\n",
    "guess_dur = guess_dur_hrs/24.\n",
    "guess_rprs = np.sqrt(guess_dep)\n",
    "\n",
    "# Limb-darkening from Claret (2017)\n",
    "guess_u1 = 0.1555\n",
    "guess_u2 = 0.4459"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit transit model\n",
    "\n",
    "# OPTIONAL: only fit data within 2 transit duration (speeds up fit):\n",
    "phase = (time - guess_epo) % guess_per\n",
    "phase[phase > 0.5*guess_per] -= guess_per\n",
    "near_transit = (abs(phase) < 2.*guess_dur)\n",
    "\n",
    "x = time[near_transit]\n",
    "y = flux[near_transit]\n",
    "yerr = flux_err[near_transit]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h2>Fit Statistics</h2><table><tr><td>fitting method</td><td>leastsq</td><td></td></tr><tr><td># function evals</td><td>114</td><td></td></tr><tr><td># data points</td><td>166</td><td></td></tr><tr><td># variables</td><td>5</td><td></td></tr><tr><td>chi-square</td><td> 185.707913</td><td></td></tr><tr><td>reduced chi-square</td><td> 1.15346530</td><td></td></tr><tr><td>Akaike info crit.</td><td> 28.6230901</td><td></td></tr><tr><td>Bayesian info crit.</td><td> 44.1830290</td><td></td></tr></table><h2>Variables</h2><table><tr><th> name </th><th> value </th><th> standard error </th><th> relative error </th><th> initial value </th><th> min </th><th> max </th><th> vary </th></tr><tr><td> per </td><td>  1.97717616 </td><td>  3.6010e-05 </td><td> (0.00%) </td><td> 1.9771744380230802 </td><td>  0.00000000 </td><td>         inf </td><td> True </td></tr><tr><td> epo </td><td>  1641.40575 </td><td>  0.00133881 </td><td> (0.00%) </td><td> 1641.4055601127118 </td><td>  0.00000000 </td><td>         inf </td><td> True </td></tr><tr><td> b </td><td>  0.99999999 </td><td>  0.20321380 </td><td> (20.32%) </td><td> 0.9 </td><td>  0.00000000 </td><td>  1.00000000 </td><td> True </td></tr><tr><td> rprs </td><td>  0.16461039 </td><td>  0.14366023 </td><td> (87.27%) </td><td> 0.1096783037143404 </td><td>  0.00000000 </td><td>  1.00000000 </td><td> True </td></tr><tr><td> ars </td><td>  6.21832487 </td><td>  0.74918675 </td><td> (12.05%) </td><td> 6.9938944820133715 </td><td>  0.00000000 </td><td>         inf </td><td> True </td></tr><tr><td> u1 </td><td>  0.15550000 </td><td>  0.00000000 </td><td> (0.00%) </td><td> 0.1555 </td><td>        -inf </td><td>         inf </td><td> False </td></tr><tr><td> u2 </td><td>  0.44590000 </td><td>  0.00000000 </td><td> (0.00%) </td><td> 0.4459 </td><td>        -inf </td><td>         inf </td><td> False </td></tr></table><h2>Correlations (unreported correlations are < 0.100)</h2><table><tr><td>b</td><td>rprs</td><td>1.0000</td></tr><tr><td>b</td><td>ars</td><td>-0.9882</td></tr><tr><td>rprs</td><td>ars</td><td>-0.9872</td></tr><tr><td>per</td><td>epo</td><td>0.9772</td></tr><tr><td>epo</td><td>ars</td><td>-0.1382</td></tr><tr><td>epo</td><td>b</td><td>0.1275</td></tr><tr><td>epo</td><td>rprs</td><td>0.1272</td></tr></table>"
      ],
      "text/plain": [
       "<lmfit.minimizer.MinimizerResult at 0x13cb0f7c0>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Fit transit model. The default is to force b to vary between [0,1], and to not fit limb darkening parameters\n",
    "best_fit = fit_transit_model(x, y, yerr, guess_epo, guess_per, guess_rprs, guess_dur, guess_u1, guess_u2)\n",
    "\n",
    "best_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Relative flux')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEGCAYAAABy53LJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA+YElEQVR4nO2deZwU1bX4v2d6pgdQFBxQESS4xpBIEMZlYpQBBHFFMQkuUYxRH3kxifqMCT/j8iQBNYmPmORFjE8NiYpbVEwwIDjD4owLCuIWFZdEFhUGFRSZ9fz+uNVDTdPd093TNdU9c76fT32m6tatW6dqquvUuefec0RVMQzDMIx0KQpbAMMwDKOwMMVhGIZhZIQpDsMwDCMjTHEYhmEYGWGKwzAMw8iI4rAF6Az69eunQ4YMCVsMwzCMguL555/fpKr948u7heIYMmQIK1asCFsMwzCMgkJE/pWo3LqqDMMwjIwwxWEYhmFkhCkOwzAMIyNMcRiGYRgZYYrDMAzDyAhTHIZhGEZGmOIwAqO2tpaZM2dSW1sbtiiGYeSQwBSHiNwhIh+KyMtJ9ouI3CIia0RktYiM8O2bIiJvessUr6yXiPxdRP4pIq+IyA1ByW50nNraWsaOHcvVV1/N2LFjTXkYRhciSIvjLmBCiv0nAAd5y8XAHwBEZA/gWuBI4AjgWhHp6x3zK1U9BDgMOFpETghGdKOjVFdX09DQQHNzMw0NDVRXV4ctkmEYOSIwxaGqS4HNKapMBOao42mgj4gMAI4HnlDVzar6EfAEMEFVt6lqldd2A/ACMCgo+XNJd+yyqaysJBqNEolEiEajVFZWhi2SYRg5IsyQIwOB93zba72yZOWtiEgf4BTgN8kaF5GLcZYMgwcPzonA2RDrsmloaCAajbJ48WIqKipCk6ezqKioYPHixVRXV1NZWdktrtkwugsFF6tKRIqBe4FbVPXtZPVU9TbgNoDy8vLQ8uMm6rLpKi/R2tralIqhoqKiy1yrYRg7CFNxrAP29W0P8srWAZVx5dW+7duAN1V1VrDi5YZYl03M4ugqXTbd1ZIyDCNcxTEPuERE5uIc4Z+o6gYRWQDM8DnExwPTAETk58DuwIVhCJwNXanLxm9hdGVLyjCM1ASmOETkXpzl0E9E1uJGSpUAqOqtwHzgRGANsA34jrdvs4hMB57zmrreKxsEXAX8E3hBRAB+p6q3B3UNuaIrdNnEWxizZs3qkpaUYRjtE5jiUNWz2tmvwPeT7LsDuCOubC0gORPQyIh4C6Ourq7LWFKGYWRGwTnHjXBI5KvpCpaUYRiZY4ojZNobmZQv7Qflqwn6+g3DyD2mOEIk6JFJuW4/1xaGjcwyChH72LEgh6ESdFiOfA/7ke/yGUY8FoPNYYojRIIOy5HvYT/yXT7DiGfOnDls374944+drhZ2SNzgpq5NeXm5rlixImwxElIoPo6gyHf5DCNGbW0tlZWVNDQ0AFBaWkpVVVW7z20hd8mKyPOqWh5fbj6OkMnEb5DNSzZTv0Rnv8htZJaRL7T37FdXV9Pc3AyAiPCd73wnrWe3K06WNcVRIHTGV0shfxkZRkdI59mPH5J+3nnnpdV2Vww7ZD6OAiETR3K2/anZOKtz1Xfb1fqAjcKiurqa+vp6mpubqa+v3+nZj1kjs2bNYvr06Rl9VMWGsmd6XD5jFkeBkO5XS0eshky/jHJloZilY4RNWVkZLS0tALS0tFBWVta6LxfPZ1frkjWLo0BI96ulI0NcM/0yytVwWhuWa+SKbC3Xuro6iorc67CoqIi6urrWffn4fIZtoZvFUUCk89XS0f7UdM4RM9vLyspy0nfbFfuAjdzjd14DOzmyO2ptl5aWJnwGg3o+Eznjkzno4689dAtdVbv8MnLkSO1O1NTU6IwZM7SmpiaQtnv27KmRSER79uyps2fPzsm5gpTZKHz8z100GtXS0tLWZzD2zMyYMUMjkYgCGolEdMaMGRmfI9kzmO3zmey4+N9RTU1NwrJEdadOndqh68wEYIUmeKeaxVFgpPNFEmR/aqIoudOmTeuQ7ND1+oC7O7ke1u1/7mK+CFVtM7w1SGs7m+czlQWUrPsr0bDd+LpA+BZ6Im3S1ZZ8sDhy8UWd7hdJkF/t2Z4rSBnNWskvgvhfp2NxxOpl8iwE+eyksoA6YnHE6nbGM08SiyP0l3pnLGErjlz9kJI9iB010TMlm4c2HRmzabczlaaRHkE9j/7nI8gPsfhzZSNfe+0nO0eqrq0wPo5McYRIrn5IyfwLs2fPDu2LJF3S+RFlowA6W2ka7VMoyjzZs5ON/KmshXz6HWZKMsVhPo4saG90Rzy5GpXhz4lRVlbGpZde2iaVa11dXaijLtrzYaTK55FtWAYbkZV/tPe/zpf4ZMmenWT+h/bCkSR6foMOKRQaibRJV1tyaXGk29caq5tL09pPoq+l2Dk6c9RFjI5+ZXbk+EL/qusOxP5HiazjdI8NapRgou6ieMu+vZGEnfH8h/GcE0ZXFS5v+IfAy0n2C3ALsAZYDYzw7ZsCvOktU3zlI4GXvGNuwYvwm2rJpeLwv7BFREVEAS0qKtLx48d3msM61cOdiULLFfGKbOrUqR3uIza6Bv5ntbi4WIuKilp/P1OnTk372M7s9vI/i/5nu6ioSEtKSnLmF4nRXpdrWPchLMVxLDAiheI4EXjcUyBHAc945XsAb3t/+3rrfb19z3p1xTv2hPbkCNriiP0QioqKWv+pufRrpDO2PN0Xd1APoL/d0tJSjUajed/HbXQO8S/e2Dqg0Wg05YdNPviwkim+XDv+U/0uZ8yY0eY901n3IZniCDwfh4gMAf6mql9JsG82UK2q93rbrwOVsUVV/8Nfz1uqVPUQr/wsf71kZJ2PY/JkeOKJeKFpbGqiqbGR4mgUgG3bttHY2EjsTu7Sqxcl0Sgfffxx62F9+/ShxKsfaye+3fjthsZG6jZtam23X79+RGNtiFDf0EB9fT2lPXoAsHHjRqeQRdhzzz0p9dWN8cmWLW3k6tO3L3123z0juZJtf15fz+fbttHY1MQnn3zSRu6yPfZI3VYm5x06FO64A3bdFSNPaWmBVatg6VLqHnmE+iVL6OHtKikpobGxsbVqz5496ek9w01NTWzZurV1X69evdi2bVvr9m69e1NcHKxrtqmpifpYzo1olOLiYpqammhsakJEApMndo6S4uKd2txeX9/mvD169EBEdqqbsI2nn4aDD85KpnzNxzEQeM+3vdYrS1W+NkH5TojIxcDFAIMHD85OutGjYa+9dmx7SrYEKPEp3M/ff59HH3mE5pYWIkVFTJwwgb332ouW999n/fr17LPPPpQkaCd++4MPPnD1Bwxgr7324pVVq1hRVweqiAjlQ4Yw/KtfdXU//JDH58+npaWFyGefccIJJxAFNmzYwIC996Z0zz13Ppcq9Rs3smzBAndcURHjjzoK+vdPKVfC7QT7egI9gU2bNrHqySdpaW5GRBi5336U7b9/eu22t93cDA89BB99BH/7G5SWYuQHbZy7d9wBt98OQNn++7NxwgTe/uQTBg4cSD3w8MMPt/5eTj/pJD4G1q1bx/sffMDbPsVx6P77079/f9a89RYHHnAAe3xlp+/PnLJhwwb++te/0uxNMow0NTHp9NMZMGBA68vyow0bWLduHT169ODt7dsZOHAgAwYM6PC5i9nxQt7gnSPW9ssrVlBbW4viulqkoQFVdffPk2/Dhg1t76tXjv/DMFckMkNyuQBDSN5V9Tfg677txUA5cAXwM1/51V5ZObDIV34MzprptK6qZHS0fz6TCUGqHRuGmoms2V7X7NmzW836nHdX3Xmn62W95prctWl0CP+zOrK0VFtEVC+8UPW995LW9w8ciR3r78YC9LTTTuvUiaMzZsxo9Vvi+WGSzTkKUq5U74Jk3WVBdOuRp8Nx1wH7+rYHeWXrcN1V/vJqr3xQgvqh09GQGYmG802bNo1Zs2bx0EMPMXz48NZhgR0Nr5CurB0JGldXV4eq0tLS0uGsZzsNUzz/fJg/H26+GS65ZIfFZISG//n9WXMz9aWl9LjxRthjj4T1/c/gzJkzW48tKioiEonQ0tJCNBpl7733zln2vFSBAmPD2cvKyigpKWkT2iPRbyvIrH7J3gXJhuLH5OvUoemJtEkuF1JbHCfR1jn+rO5wjr+Dc4z39db30MTO8RPbkyHsCYDpkOorI5HzPXZMkKOQOmrV5OKLLGk7r72mWlSkevnlWbVr5JbY/+nIoiJV0H9feGHGxyYa7hrUc+Qfsh4/Umr27Nk6depUnTp1atLzdbbFkahObHhz/Gz1XL4TCMPiEJF7cZZDPxFZC1yLcxGgqrcC83Ejq9YA24DvePs2i8h04DmvqetVdbO3/p/AXbju9Me9pdPJ9WSdRJOmYl9i/gQz2U4uyoZ0vmCS3Yf4yYp+aykT4r++5syZs+N8U6bQ8rvf8b8lJYycODH/J011YWL/772mTKHxww/Z99e/zvjYZL+nVPvSJVWgQBFpDZ4YC9z5hz/8oUMyd4R02o6VpQqi6K+XcxJpk6625Nri6Kwx1e1ZHNm2masgcOl+GeVqYlT8/JS7Z8zQBtCbRWzYbz6wdKnzPf3yl4GdItsv6mQWfbaTEjuD9q412SRg/+8lldWUDlisqtzRmWPLk5mk2baVyx9IOvchF/cqdg/iuxcOPPBAfQB0E2gPaHcymREgLS2qxxyjuvfeqp99FsgpcvERkut8G0GR6UdZaWlpa/eaf4CBdPCjyhRHDglrFmdHybXCS9U3naxOLhSff+bxGDdYV89OMpnM6CSeeca9Tm65JavD03lx58NkwEQEoXTSvdaamhqdOnVq64TbmFXuHxnWkXtliiPH5NMXSrqyBKHw0jH3c3mvampqdPz48Tu670DfAK3Os5dJt+Oyy1SjUdWPPsr40HSfy3z8YOuMSAzttZsoasTUqVNThh1KF1McOSB+7Hk+KI5MH9yg5O7s7ju/iX5NJKIKemBpaej/j25Jc7PqoEGqp5ySslqyZy+TZ6cz5iBlQqbPfSYyzZ49W8ePH6+zZ89ut81oNKoi0sbqzsX1m+LoIKmctGG+rPIliVMmCiwXD7S/jRfuu08V9J0f/CDr9owOsHy5e5X85S9Jq6R6PoKYLNrZA1jSfe6DqhvrnirN8cdTMsUR9gTAvMY/1DSdnMeZtpmLoXKdOekn1YTAdIcnpjOpMJ171GYockUF3HQTQ55+OncXa6TPffe50C+nnpq0SrIJc7W1tVx66aUuBE4kwqxZs3KeozzRbzRXv8NMhuWmM2kwJte///3vtCcYVldX09TUhKrS1NSU08mISUmkTbrako3Fkcjx21GLI8j+0M7oNsuFdRNY+OibbnJfvWvWZCyT0QGamtxIqtNPT1ktvjslNpAiqNwxqZ6jsHwl7Z03216NIK8HszgyI/7rYOXKlUyZMgWA8847r7VOJl8sQYUpCHoiYIxcWDfttZH1PfrWt+DKK+H++2HatIzlMrJk+XJ4/30XSbodxItwrKr88Ic/pKmpieLiYiKRCJA8vEc2pLIEggwXkq1M8XIBXHTRRQwePLjdd0yQkxGTkkibdLWloxZHrnwa+TgqJFNy7Z9ItC/re1RRoTpsWNZyGVnwve+p9uyp+umnKaslS4CWKndMUOTr7zAf5SKsfBz5QLb5OPz9jX/84x9pbm4mEokwffp0piX5qm2v77Sg8gqHRNb36H/+By6/HN56C/xh3I1gaGqCffaBykpn6aXA79uKRCKICE1NTRkHz8wV+fo7zDe5kuXjMMWRBulGie1INFmj47zwwAOM+Na3eOdHP2K/WbPCFqfrs2gRjBsHDz4IZ5zRbvX46LTpvCDz7UXa3cjXRE4FQbp9iGH1nRqe0p4yhWeBTb/9Le9Pnmz3PmgeeAB22QVOPDGt6vG+uPb+P/Yhlr8UhS1AoVBRUcG0adNSPrgxx28kEgk+Hr7RhpjSfhT4eksLT8+fH7ZIXRtV6h95hH8OGULtqlWBnCLRh5iRH5jiyCExy2T69On2ddTJxJT230UoBg59//2wRerSrJo7l9IPP2TWq68yduxYamtrc34O+xDLX0xx5Jh0LBMj91RUVDBr1iyeLypiM7D+zjsDeZkZjo133w3AP1QDswbsQyx/MR+H0WWoq6ujGVgEjG1uZk5Vlb1sAqJ882beEGFtUVGg1kBnzVEyMsMUh9FliHVtPLF9O99S5YTBg8MWqWtSX0/fF19k+xlnMH3ECBvx1A0xxWF0GWJdG88/+ijceCPDP/ggbJG6Jk89Bdu2MeC885h2yilhS2OEQKA+DhGZICKvi8gaEflpgv1fEJHFIrJaRKpFZJBv340i8rK3TPaVjxWRF0RklYgsF5EDg7wGo7CoqKjgkhtugKFDYcGCsMXpmixcCMXFbuKf0S0JTHGISAT4PXACMBQ4S0SGxlX7FTBHVYcB1wMzvWNPAkYAw4EjgStEZDfvmD8A56jqcOAe4GdBXYNRwIwfD0uXwrZtYUvS9Vi4EI4+Gnr3DlsSIySCtDiOANao6tuq2gDMBSbG1RkKPOmtV/n2DwWWqmqTqn4GrAYmePsUiCmR3YH1AclvFDLHHw/19U55GLnjww9h5UqnmI1uS5CKYyDwnm97rVfm50Vgkrd+OtBbRMq88gki0ktE+gGjgX29ehcC80VkLXAucENA8huFzLHHQmkp6++6i5kzZ9rQ3FyxaJH7a4qjWxP2PI4rgFEishIYBawDmlV1ITAfqAHuBWqBZu+Yy4ATVXUQcCdwc6KGReRiEVkhIis2btwY8GUYeUevXnw8bBif3H8/V199dWCT1LodCxZAWRkcdljYkhghEqTiWMcOKwFgkFfWiqquV9VJqnoYcJVX9rH39xeqOlxVxwECvCEi/YGvquozXhP3AV9LdHJVvU1Vy1W1vH///rm8LqNAeL5fP76kygALWZEbVJ1/47jjwMuhYXRPglQczwEHich+IhIFzgTm+SuISD8RickwDbjDK494XVaIyDBgGLAQ+AjYXUQO9o4ZB7wW4DUYBUzZ2WcDcIKIhazIBS+/7JI2HX982JIYIRPYPA5VbRKRS4AFQAS4Q1VfEZHrcclB5gGVwEwvqctS4Pve4SXAMi9j2Bbg26raBCAiFwEPiUgLTpFcENQ1GIXN8HPO4dMf/IBvl5ZSfv31Nkmtoyxc6P6OGxeuHEboWD4Oo8tSW1vLW8ccw/HNzQzp0YNFTz5pyqMjHH88rF0Lr7wStiRGJ5EsH0fYznHDCIzq6mqebGmhP3Cw+Tg6xuefu6HNNprKwBSH0YWprKzkqWgUgLGRiPk4OsLy5bB9u/k3DMAUh9GFqaio4K6qKj7q25f/d9RR1k3VERYsgGjUzY8xuj2mOIwuTUVFBX0nTWKPl16C5ub2DzASs2gRfP3r0KtX2JIYeYApDqPrM3o0fPwxBJTitCtRW1u780z7zZth9Wp3Hw0DC6tudAdiL7yqKhg5MlxZ8pja2lrGjh1LQ0MD0Wh0R9a9Zcvc5D/zERkeZnEYXZ999oGDD3aKw0hKdXU1DQ0NNMfPtK+uhh494PDDwxTPyCPaVRwi8t247YiIXBucSIYRAGPGuOGkjY1hS5K3xDIoRiKRtjPtlyyBr30NSktDlc/IH9KxOMaKyHwRGSAiXwaeBiwQv5GXJOyjB9dd9emn8Pzz4QhWAMQyKE6fPr21m+rZBQvQVat4b//9wxbPyCdUtd0FmAxsAv4FHJ3OMfm0jBw5Uo2uT01Njfbs2VMjkYj27NlTa2pqduz84ANVUJ0xIzwBC4yamho9IxpVBT0uGm17P41uAS481E7v1HS6qg4CfgQ85CmOc0XExuQZeUfSPnqAPfeEr3zF/BwZUF1dzdGNjWwHnmps5LrrrrPQ9AaQXlfVY8A1qvofuJwZb+Ii3xpGXpG0jz7G6NFuBnRDQyjyFRqVlZWMwiXD+VyVRYsWWV4TA0hPcRyhqosAPOvl17hsfYaRVyTqo29DZaWLufScffekQ8WXvsRhIqzdf3+KiopoaWmxvCYGkN48jtO88ObxvJFjWQyjw1RUVCQPLRILl7FkCRx9dOcJVagsX460tDDi8ssp/fGPW+d3WMwvIx3F4R+83QMYC7wAzAlEIsMIiNo33+SAvfcm+sgjvDZ6NNXV1VRWVloMq2QsWQKlpXz5u99l8YgRdr+MVjLOxyEifYC5qjohEIkCwPJxGLFZ0Tdt3853VNmrpITtLS1tZ0gbbTn8cBebasmSsCUxQiKX+Tg+A/bruEiG0XnERlxVqbIL8NXGxsSjrwzHli3wwgsWZsRISLtdVSLyGBAzS4qAocD9QQplGLkmNuLqqfp6aGlhdCTCM2B99slYvhxaWmDUqLAlMfKQdHwcv/KtNwH/UtW1AcljGIExZcoUALYtXMgV/fqxy2mnWZ99MpYscfk3jjoqbEmMPKRdxaGqWXdwisgE4DdABLhdVW+I2/8F4A6gP7AZ+HZMKYnIjcBJXtXpqnqfVy7Az4FvAs3AH1T1lmxlNLo+8VFfrz3xRPZ+/HGmLV8OJSVhixc6tbW1rY5vcN16l/ztb/Q+4gjLv2EkJKniEJGtuC4qYUdXFbFtVd0tVcMiEgF+D4wD1gLPicg8VX3VV+1XwBxV/ZOIjAFm4mamnwSMAIYDpUC1iDyuqluA84F9gUNUtUVE9szkgo3uR/yM8ppolEnbtsGKFdDNrQ2/Uo1EIogIPRob+XFLC2vPP59BYQto5CWpnONfVdXdVLW393c3/3YabR8BrFHVt1W1AZgLTIyrMxR40luv8u0fCixV1SZV/QxYDcRGcX0PuF5VWwBU9cM0ZDG6MfEzyvc95xwA/jVnTuKAiN0Iv1JtbGykoaGBo1paKAaWFVnWBSMxqZ6MBwBEZHGWbQ8E3vNtr/XK/LwITPLWTwd6i0iZVz5BRHqJSD9gNM7KADgAmCwiK0TkcS+W1k6IyMVenRUbN27M8hKMrkD8jPLDTzqJbfvtxxu33cbVV1/drcNo+JVqSUkJ0WiU0SI0APt7CtYw4knl4ygSkf8HHCwil8fvVNWbc3D+K4Dficj5wFJgHdCsqgtF5HCgBtiIC5cTSxhdCmxX1XIRmYTzkRyTQL7bgNvAzePIgaxGARM/o/y1Pfek4p13EGgdktsdneQxper3cQyePJntffpw5Jgx4Qpn5C2pFMeZwGlenWzyb6xjh5UAMMgra0VV1+NZHCKyK3CGqn7s7fsF8Atv3z3sCHGyFvirt/4wcGcWshndnF1POYVdn3mGI4qKWNnNh+S2Uaqffgrr18O554YrlJHXJFUcqvo6cKOIrFbVx7No+zngIBHZD6cwzgTO9lfwuqE2e/6KaTjrIeZY76OqdSIyDBgGLPQOewTXdfUOLlqvxcwyMuaLF10EP/sZPx83jh7XXtstrY2E1NRAc7NN/DNSks5w3GyUBqraJCKXAAtww3HvUNVXROR6XHKQeUAlMFNEFNdV9X3v8BJgmRdccQtumG6Tt+8G4G4RuQz4FLgwG/mMbs6ee8LQoYwW6fYjq9pQXQ3FxS5VrGEkIZ0JgFmjqvOB+XFl1/jWHwQeTHDcdtzIqkRtfsyO+R2GkT2jRsGf/+zykHfj+Rz+eRwVS5a4GFW77BK2WEYeE6jiMIy8prIS/vAHF5PpyCPDliYU/PM4+pSU8GFjI0U//nHYYhl5TjqpY3uJyNUi8kdv+yAROTl40QwjYGJxmLpx9Ff/PI6RDQ0UmX/DSIN0ZvjcCdQDsY7gdbiQH4ZR2Oy1F3zpS65fv5vin8cxpqgIjUTMv2G0SzpdVQeo6mQROQtAVbdJkpSAhlFwjBoFf/kLNDU5p3A3wz+P43tz5yI9ekDvbEbfG92JdCyOBhHpiRevSkQOwFkghlH4VFa6uQsrV4YtSWhUVFQw7Yc/ZLfXXoPRo8MWxygA0lEc1wH/APYVkbuBxcCVQQplGJ1GzM/RjburAHjqKTe6zBSHkQbtKg5VXYib3X0+cC9QrqrVwYplGJ3E3nvDIYeY4qiqcl11Rx8dtiRGAZDOqKrHgPFAtar+TVU3BS+WYXQio0bBsmXOz9FdqaqCI46AXXcNWxKjAEinq+pXuCCCr4rIgyLyDRHpEbBchtF5VFbC1q2walXYkoTD1q0uN0lcN1VtbW23DztvJCbdDIBLvPhRY4CLcDGl0snJYRj5j9/PUV4eqiihsGwZNDfzyl57MW/mzNaAj/6siYsXL7Z4XkYraY0/9EZVnQJMxmXm+1OQQhlGZ1L77rsc3K8fRQ8/TN8rrghbnM6nqoqWkhKOufJKtjQ2Eo1GmTJlSpusid017LyRmHR8HPcDr+Gsjd/h5nX8IGjBDKMziIXceKiujqKaGmqXLw9bpM6nqoq1AweypbGxVVEAbbImduew88bOpOPj+D+cspiqqlWxlK2G0RWIhdyoUmV34J9z54YtUufy8cewciUyZkwbRXHeeee1yZpo1obhJ2lXlYiMUdUngV2AifGTxVX1rwkPNIwCIhZyY3l9PbS0uDDr3YmlS6GlhX3PO4/FF164I0qupyhMYRiJSOXjGAU8ifNtxKPsyMJnGAWLP+TG57feypB33wXiQo135ZdnVRX06AFHHUVFaWnXvlYjZ6TKAHitt3q9qr7j3+dl9TOMLkFr6tR33oH77qN2+XLGjh/fPUYUVVW5oIalpWFLYhQQ6fg4HkpQtlPyJcMoeCorYcsWXps7d6cRRV1yTkNdHbz4ooUZMTImlY/jEODLwO4iMsm3azfAJgAaXQ9vPseYoiKi0WirxVFWVlbQcxqSdrvF8pCY4jAyJJWP44vAyUAf2vo5tuImARpG12LgQDjoIIa8+26r36OysrJNsqNCm9Pgz/AXjUaZNWsWdXV1TolUVUGvXi5VrGFkQCofx6PAoyJSoapZ2eciMgH4DRABblfVG+L2fwE3C70/sBn4tqqu9fbdyI7c4tNV9b64Y28BLlBVC65j5I5Ro+CBB6h4+OE2ysFvgRTSnAa/0quvr+eSSy6hpaWFaDTKpr33ptfXvw7RaNhiGgVGOjPHV4rI93HdVq1dVKp6QaqDvBAlvwfGAWuB50Rknqq+6qv2K2COqv5JRMYAM4FzReQk3Az14UApUC0ij6vqFq/tcqBvmtdoGOlTWQm33+76/keMANqOvCq0UVax4cYNDQ2ICM3NzbS0tNCnvp5e77wDF18ctohGAZKOc/zPwN7A8cASYBCuu6o9jgDWqOrbqtoAzAUmxtUZihvyC1Dl2z8UWKqqTar6GbAamACtCumXWE4QIwjGjHF/Fy1qU1xRUcG0adMKSmnADqU3ffp0fv/731NaWkokEmFsJOIqmH/DyIJ0LI4DVfWbIjLRswzuAZalcdxA4D3f9lrgyLg6L+JyffwGOB3oLSJlXvm1IvJroBcwGohZKpcA81R1Q6oMtiJyMXAxwODBg9MQ1zCAAQPgy192iuPKrvFt0jrcGDj00EOprq7mO88+C4sXw8iRIUtnFCLpKI5G7+/HIvIV4H1gzxyd/wrgdyJyPrAUWAc0q+pCETkcqAE2ArVAs4jsA3wTqGyvYVW9DbgNoLy8XHMkr9EdGDcObr0Vtm93k+O6EK1K5ItfhGOO6ZZ51o2Ok05X1W0i0he4GpiH+/K/KY3j1gH7+rYHeWWtqOp6VZ2kqocBV3llH3t/f6Gqw1V1HCDAG8BhwIHAGhF5F+glImvSkMUw0ue445zSeOqpsCUJhvXr4Y03rJvKyJp08nHc7q0uAfbPoO3ngIO8WebrgDOBs/0VRKQfsNkLnDgNN8Iq5sfoo6p1IjIMGAYsVNUmnL8ldvynqnpgBjIZRvsce6z7El+0CMaODVua3POk51Y0xWFkSaoJgJenOlBVb25nf5OIXAIswA3HvUNVXxGR64EVqjoP1+U0U0QU11X1fe/wEmCZ58PYghum243zehqdSu/eUFEBTzwBM2eGLU3ueeIJKCuD4cPDlsQoUFJZHL072riqzgfmx5Vd41t/kAThS1R1O25kVXvt2xwOI6fEZll/+4tfZN//+z8XlqOsLGyxcocqLFzo/DixkVWGkSGpJgD+d2cKYhhh459lvSASoVqV12+9lb8WFRXc/I2kvPQSvP8+jB8ftiRGAdOuj0NEDgb+AOylql/xfA6nqurPA5fOMDoR/yzrWlU+j0ZZfu21XA0FGaMqIQsXur/jxoUrh1HQpDOq6o84x3UjgKquxjm6DaNLEZtlHYlEiJSW8vo++zC6ublNjKqCZ8ECN09l0KCwJTEKmHQURy9VfTauzBzVRpfDP8t68eLF7D5pEvsDB3rRcgslRlXSEPDbtsGyZdZNZXSYdGb/bBKRA3BZ/xCRbwAbApXKMELCP8uaPfaAm2/mllNOoc9PflIQ3VTx0XDbdK8tWwb19aY4jA6TjsXxfWA2cIiIrAMuBaYGKZRh5AUHHwz77ssJkUhBKA0gYQj4VhYscJn+jj02NPmMrkG7isMLUngcLvT5Ibhc5F8PWjDDCB0ROP54NxGwsbH9+nmA30+zU/fawoVwzDHUvvhi18tmaHQqqSYA7oazNgYCjwKLvO3/wkWrvbszBDSMUJkwwYVZf/ppF9spz0kaAn7dOnjlFf5VWVnQ2QyN/CCVj+PPwEe4AIMX4WJJCXC6qq4KXjTDyAOOO85NlHv88VbFkTQVa57Qxk8TwxuG+0RRUcFmMzTyh1SKY39VPRRARG7HOcQHe7O6DaN7sPvu8LWvOcUxY0Zq53Mes+mee+jRuzctX/5ywWYzNPKHVD6O1k5dVW0G1prSMLolJ5wAq1bBhg2pnc95ytPLlhFZtIgHP/2USy+7jFmzZrUOOS4EpWfkH6kUx1dFZIu3bAWGxdZFZEtnCWgYoXPCCe7v44+ndj7nKW/NmUNfYJ4qDQ0N1NXVFWQ2QyN/SBWryiKgGQbAV7/qZlo/9hgVF1xQcPnHR2/dSj2w2JvIWFZWxsyZMwtGfiP/sPRfhtEeInDqqXDXXfD554mdz/mKKvs8/zwfHXUUPz31VMrKyrj00ksLzkdj5BfpTAA0DGPiRBeyY9GisCXJjNdfhzVr6HvuuUybNo26urqC89EY+YcpDsNIh8pKl+Bp3rywJcmMmLwnnwy0M0HQMNLEuqoMIx2iUeckf+wxaGmBogL55nrsMZfpb/BgIMUEQcPIAFMchpEuEyfC/ffDM8+41LL5zqZNUFMDV13VprigfDRGXlIgn02GkQeceCIUF8Ojj4YtSXrMn++so1NOCVsSo4sRqOIQkQki8rqIrBGRnybY/wURWSwiq0WkWkQG+fbdKCIve8tkX/ndXpsvi8gdIlIS5DUYRit9+sCoUYXj53jsMRgwAEaOTJ6jwzCyIDDFISIR4PfACcBQ4CwRGRpX7VfAHFUdBlwPzPSOPQkYAQwHjgSu8IIugguueAhwKNATuDCoazCMnZg4EV57Dd58M2xJUtPQ4MKon3wytc88w9ixY7n66qsZO3asKQ+jwwRpcRwBrPHCsjcAc4GJcXWGAk9661W+/UOBparapKqf4aLxTgBQ1fnqATwLWA5Mo/M49VT3N9+7qxYvhq1b4dRTCzJMipHfBKk4BgLv+bbXemV+XgQmeeunA71FpMwrnyAivUSkHzAa2Nd/oNdFdS7wj0QnF5GLRWSFiKzYuHFjhy/GMAD4whfcTPJ8VxwPPgi77QbjxtkQXCPnhO0cvwIYJSIrcQmi1gHNqroQmA/UAPfiQrs3xx37vzirZFmihlX1NlUtV9Xy/v37B3YBRjdk0iR46imX4yIfaWyERx5x1lFp6U651G1EldFRglQc62hrJQzyylpR1fWqOklVD8Pl+0BVP/b+/kJVh6vqOFwekDdix4nItbiMhJcHKL9hJOass0AV5s4FyD/Hc1UVbN4M3/hGa1FFRYUFNjRyRpDzOJ4DDhKR/XAK40zgbH8Frxtqs6q2ANOAO7zyCNBHVetEZBgwDFjo7bsQOB4Y6x1nGJ3LQQfB4YfDPfdQ+7Wv5V9+jvvvh113hfHjw5XD6LIEZnGoahNwCbAAeA24X1VfEZHrRcTzMFIJvC4ibwB7Ab/wykuAZSLyKnAb8G2vPYBbvbq1IrJKRK4J6hoMIylnnw0vvMDq++/PL8fz55/DAw/AGWdAz57hymJ0WQKdOa6q83G+Cn/ZNb71B4EHExy3HTeyKlGbNtvdCJ/Jk+G//ouTPvmEy/Ipo968ebBlC5x3XrhyGF0aewkbRjYMGABjxjBo6VIWL1pE9ZIl+RH7ac4c2HdfF5TRMALCFIdhZMvZZ8MFF1BRXEzFtGlhSwPvv+8m/V15ZeEEYTQKEnu6DCNbJk2C0lK4++6wJXHcey80N8O554YtidHFMcVhGNmy++4uz8V990GTG7sR6tDcOXPcaK8vfanzz210K6yryjA6wjnnwEMPwZNPUtu7d3hDc1evhlWr4JZbWotqa2st74YRCKY4DKMjnHAC7LEH/N//UT18+E5Dczvthf3nP7uQ72eeCTilkXfzS4wug3VVGUZH6NHDDX19+GGOGzYsnJhQTU3Oz3LiieCF17HAhkaQmOIwjI5y8cXQ2Mjhr74aaEyopP6TRx+FDRvg/PNbiyywoREk4qKTd23Ky8t1xYoVYYthdGWOPRbWr4fXX4dIJOfNp+x6OvpopzjefLPNuc3HYXQUEXleVcvjy83iMIxc8MMfwltvBRZuPWnX09NPQ00N70ycyMybbmpjjVhgQyMozDluGLng9NPhgAPgppvcukhOm491Pe0U2uTmm2nadVeOuPVWPmpsNEe40SmYxWEYuSASgcsvh2eegeXLc958wpwa774LDz3EcyNG8FFjoznCjU7DLA7DyAG1tbUs37iRS/v0Yeu0acw+6aSc+xYqKiratnfLLVBURPTyy4k+91z+BFo0ujzmHDeMDuJ3XF8LXN3czKFFRbxVWrpTt1G8wzprB/Ynn7hghqeeCn/5iznCjUBI5hw3i8MwOojfcX0L8F/AZS0tXBw3CTB+ZNSsWbO49NJLs5ukd/vtsHUrXHYZkMAaMYwAMcVhGB3E77jeWlzMnxob+W5LCzcVF1NZWdlqDfz73/9uMzLqoYceym6m+ZYtcMMNMHYsjBwZ/AUaRhymOAyjg8Qc17GuougHHyDf/CbLjzmGN6HVyohEIhQXu59cNBrljDPOYNmyZZn7Jn79a9i0CWbODOyaDCMVpjgMIwfs1FU0bRr9pk9n/gEHtFoVABdddBGDBw9u9UUceuihmfk83n/fKY5vftNFwjWMEDDnuGEEwaefwhe/yKe7785e77xDfRpzLNIKTHjuuS6M+yuvwEEHdcKFGN2ZUGaOi8gEEXldRNaIyE8T7P+CiCwWkdUiUi0ig3z7bhSRl71lsq98PxF5xmvzPhGJBnkNhpEVu+4KM2ey62uvceuRR3LRRRe16/xuNzDhE0/AX/4CP/mJKQ0jVAJTHCISAX4PnAAMBc4SkaFx1X4FzFHVYcD1wEzv2JOAEcBw4EjgChHZzTvmRuB/VPVA4CPgu0Fdg2F0hNoDDmBZURGnLFnCorvuard+ysCEW7fC977nFMZVVwUms2GkQ5AWxxHAGlV9W1UbgLnAxLg6Q4EnvfUq3/6hwFJVbVLVz4DVwAQREWAM8KBX70/AacFdgmFkT/XSpXwXKAH+tH07y594ImX9hLPDAbZtc5kG330X/vhHF8rdMEIkSMUxEHjPt73WK/PzIjDJWz8d6C0iZV75BBHpJSL9gNHAvkAZ8LGqNqVoEwARuVhEVojIio0bN+bkggwjEyorK1lbWsp3ioo4CrigqsrlBE/BToEJ6+td7Ktly1w31ahRwQtuGO0QdqyqK4BRIrISGAWsA5pVdSEwH6gB7gVqgdS/uDhU9TZVLVfV8v5echvD6ExiFsRhP/85//rRjyirrnYT9tIdkNLYCJMnw8KFbsKfl93PMMImyOG463BWQoxBXlkrqroez+IQkV2BM1T1Y2/fL4BfePvuAd4A6oA+IlLsWR07tWkY+USbYbrFxW4o7aBBcOWVqQ9sbnaJmR591MWkuuCCwGU1jHQJ0uJ4DjjIGwUVBc4E5vkriEg/EYnJMA24wyuPeF1WiMgwYBiwUN3Y4SrgG94xU4BgEiAYRgdImK3vppucBfGTn7hhtZs3Jz64qQmmToV77nGT/H7wg/TaN4zOQlUDW4ATcZbCW8BVXtn1wKne+jeAN706twOlXnkP4FVveRoY7mtzf+BZYA3wQOyYVMvIkSPVMDqLmpoa7dmzp0YiEe3Zs6fW1NTs2NnQoHrNNaqRiOqAAao33KD6wguq69ervv666u9+pzpkiCqoXnVV5u0bRg4BVmiCd2qgM8dVdT7OV+Evu8a3/iA7Rkj562zHjaxK1ObbuBFbhpGXJJqP0dpdVVIC//3fMHEiXHop/PSnbvFz5JGue+rkkzNv3zA6AQs5Yhg5Jmm2Pj8jRsDSpbBmDaxe7UKJ9OgB5eVw6KEpMwim1b5hBIgpDsPIMfFBD1NaAwceCAceuCNG1WefUdFO2tmM2jeMALBYVYYRMmnFqDKMEAglVpVhGO3Tbowqw8gzrKvKMEIi1j1VVlZmPgujoDDFYRghkCiNbF1dnfksjILAFIdhhEB891RdXR3Tpk0LWyzDSAvzcRhGCKQMoW4YeY5ZHIYRAtkOqW03taxhdAKmOAwjJHbKU94ONmzXyBesq8owCgQbtmvkC6Y4DKNAML+IkS9YV5VhFAgWasTIF0xxGEYBkalfxDCCwLqqDMMwjIwwxWEYhmFkhCkOwzAMIyNMcRiGYRgZYYrDMAzDyIhAFYeITBCR10VkjYj8NMH+L4jIYhFZLSLVIjLIt+8mEXlFRF4TkVtEXFo0ETlLRF7yjvmHiPQL8hoMIwxqa2uZOXMmtbW1YYtiGDsR2HBcEYkAvwfGAWuB50Rknqq+6qv2K2COqv5JRMYAM4FzReRrwNHAMK/ecmCUiCwHfgMMVdVNInITcAlwXVDXYRidjYUWMfKdIC2OI4A1qvq2qjYAc4GJcXWGAk9661W+/Qr0AKJAKVACfACIt+ziWSC7AesDvAbDCJREloWFFjHynSAnAA4E3vNtrwWOjKvzIjAJZ0WcDvQWkTJVrRWRKmADTlH8TlVfAxCR7wEvAZ8BbwLfT3RyEbkYuBhg8ODBubomw8gZySyLWGgRywho5CthO8evwHVBrQRGAeuAZhE5EPgSMAingMaIyDEiUgJ8DzgM2AdYDSTMfqOqt6lquaqW9+/fvxMuxTAyI5llEQstMn36dOumMvKSIC2OdcC+vu1BXlkrqroeZ3EgIrsCZ6jqxyJyEfC0qn7q7XscqAC2e8e95ZXfD+zkdDeMQiCVZWGhRYx8JkiL4zngIBHZT0SiwJnAPH8FEeknIjEZpgF3eOv/xlkixZ6VMQp4Dad4hopIzIQY55UbRsFhloVRqARmcahqk4hcAiwAIsAdqvqKiFwPrFDVeUAlMFNEFFjKDn/Fg8AYnC9DgX+o6mMAIvLfwFIRaQT+BZwf1DUYRtCYZWEUIqKqYcsQOOXl5bpixYqwxTAMwygoROR5VS2PLw/bOW4YhmEUGKY4DMMwjIwwxWEYhmFkhCkOwzAMIyNMcRiGYRgZ0S1GVYnIRtzQ3XyjH7ApbCEywOQNlkKSt5BkBZM3W76gqjuF3ugWiiNfEZEViYa65Ssmb7AUkryFJCuYvLnGuqoMwzCMjDDFYRiGYWSEKY5wuS1sATLE5A2WQpK3kGQFkzenmI/DMAzDyAizOAzDMIyMMMVhGIZhZIQpjgAQkT1E5AkRedP72zdJvSlenTdFZIpX1ltEVvmWTSIyy9t3vohs9O27MGx5vfJqEXndJ9eeXnmpiNwnImtE5BkRGRKmrCLSS0T+LiL/FJFXROQGX/2c3lsRmeDdkzUislOysVT3RkSmeeWvi8jx6bYZhrwiMk5EnheRl7y/Y3zHJHwuQpZ3iIh87pPpVt8xI73rWCMit4iIhCzrOXHvghYRGe7tC+zepoWq2pLjBbgJ+Km3/lPgxgR19gDe9v729db7Jqj3PHCst34+Lv96XskLVAPlCY75T+BWb/1M4L4wZQV6AaO9OlFgGXBCru8tLv/MW8D+3nleBIamc2+AoV79UmA/r51IOm2GJO9hwD7e+leAdb5jEj4XIcs7BHg5SbvPAkcBAjweezbCkjWuzqHAW0Hf23QXsziCYSLwJ2/9T8BpCeocDzyhqptV9SPgCWCCv4KIHAzsiXvBBUlO5G2n3QeBsTn4istaVlXdpqpVAKraALyAS2mca44A1qjq29555npy+0l2byYCc1W1XlXfAdZ47aXTZqfLq6or1aWABngF6CkipTmSK+fyJmtQRAYAu6nq0+rezHNI/GyFJetZ3rF5gSmOYNhLVTd46+8DeyWoMxB4z7e91ivzE/v68A99O0NEVovIgyKyL7khF/Le6ZnMV/se+tZjVLUJ+AQoywNZEZE+wCnAYl9xru5tOv/bZPcm2bHptBmGvH7OAF5Q1XpfWaLnImx59xORlSKyRESO8dVf206bYcgaYzJwb1xZEPc2LQJLHdvVEZFFwN4Jdl3l31BVFZcaNxvOBM71bT8G3Kuq9SLyH7ivlDEJj+xcec9R1XUi0ht4yJN5ToZttBL0vRWRYtyP8BZVfdsrzvreGiAiXwZuBMb7inP6XOSIDcBgVa0TkZHAI57seYuIHAlsU9WXfcWh3ltTHFmiqscl2yciH4jIAFXd4JnAHyaotg6Xcz3GIFy/ZayNrwLFqvq875x1vvq34/r7Q5dXVdd5f7eKyD0483yOd8y+wFrvZb074L+GTpfV4zbgTVWd5Ttn1vc2yfn9FssgryxRnfh7k+rY9toMQ15EZBDwMHCeqr4VOyDFcxGavJ71Xu/J9byIvAUc7NX3d1vm6v526N56nEmctRHgvU0L66oKhnlAbNTRFODRBHUWAONFpK+4kUHjvbIYZxH3sHgvyhinAq+FLa+IFItIP0++EuBkIPZl5G/3G8CTcd1unSqrJ+PPcT/MS/0H5PjePgccJCL7iUgU98Ofl+I6/PdmHnCmN9JmP+AgnNM2nTY7XV6vy+/vuAELT8Uqt/NchClvfxGJeHLtj7u/b3vdn1tE5Civ2+c8Ej9bnSarJ2MR8C18/o2A7216hOWV78oLrn9yMfAmsAjYwysvB2731bsA5/xcA3wnro23gUPiymbiHJAvAlXx+8OQF9gFN/JrtSfbb4CIt68H8IBX/1lg/5BlHQQoTims8pYLg7i3wInAG7gRNVd5ZdcDp7Z3b3Bdcm8Br+Mb2ZOozRw+s1nJC/wM+Mx3P1fhBnQkfS5ClvcMT55VuMERp/jaLMe9gN8CfocXWSMsWb19lcDTce0Fem/TWSzkiGEYhpER1lVlGIZhZIQpDsMwDCMjTHEYhmEYGWGKwzAMw8gIUxyGYRhGRpjiMDJCRMp8ETnfF5F13vqnIvK/AZyvv7iIoSt94SHyFhEZLiInBtj+vV5YlMviyk8TkaFBndd3nttj5xGR/5einojIkyKyW4J914nIFTmU6WQRuT5X7RntY4rDyAhVrVPV4ao6HLgV+B9ve1dV/c8ATjkWeElVD1PVNsEeYxO58ozhuHH7O+HNCs4aEdkbOFxVh6nq/8TtPg0XWTfn5/Wjqheq6qveZlLFgbsHL6rqllydOwV/B04RkV6dcC4DUxxGjhCRShH5m7d+nYj8SUSWici/RGSSiNwkLtfBP7zZrrH8B0vE5XFYEDd7G3G5B24CJnpWTU/Psvm1iLwIVIjI5SLysrdc6h03RFzOjbtE5A0RuVtEjhORp8Tl5zgigfwREfmV185qEflBKhnF5UO4UUSe9c5xjDcz+HpgsifvZO9e/FlEngL+7FlQD4nIc95ydAJZeojInd79Wikio71dC4GBXtvH+Op/DTfb/ZfevgM8+WaJyArgRyJyis9yWyQie/n+V3d49d8WkR965buIy13yondPJvuuu1xcLpOe3vnuTvBInINv5rWIXOXdp+XAF33lF3n34UXvvvQSl5PmHd9zsltsW0R+KCKvev+jueBiluFCypycQA4jCDpztqEtXWsBrgOu8NYrgb/5ypcDJcBXgW3syHvxMO7ruASoAfp75ZOBOxKc43x8eTJwM7+/5a2PBF7CzaTdFTeL9jBczoUmXA6DItws2ztweRYmAo8kOM/3cCGti73tPVLJiHtR/dpbPxFYlETe67zz9/S27wG+7q0PBl5LIMt/+c5zCPBv3OziISTPJXEX8A3fdjXwv77tvtA64fdCn+zXeddYCvTDxUgqwc2w/qPv+N197ZZ765+meDb+BfSO+z/1AnbDzZCOPTdlvmN+DvzAW78TOM1bv9gn73qg1Fvv4zv2HOC3Yf8mustiQQ6NoHhcVRtF5CVcMpt/eOUv4V6AX8Ql/nlCXEToCC5yaXs046KBAnwdeFhVPwMQkb8Cx+Bi/7yjqi955a8Ai1VVPXmGJGj3OFwynSYAVd0sIl9pR8a/en+fT9JmjHmq+rnvPENlRxTs3URkV1X91Ff/68BvPTn+KSL/wgXiy7Tb5z7f+iDgPs9iigLv+Pb9XV0o9HoR+RAXqv4l4NciciPugyDTnDB7qOpWb/0Y3P9pG4CI+GM1fUVc/LA+OOUfi9d2O3Al8AjwHeAir3w1cLeIPOLti/EhsE+GMhpZYorDCIpYBNIWEWlU77MQaME9dwK8oqoVGba7XVWb0z2/75z1vvV0n/v2ZIy12dxOm5/51ouAo1R1e5oydAT/eX8L3Kyq80SkEmdpxPDfq2ac1fWGiIzAWVM/F5HFqpqJA7pJRIpUtaWdenfhLIsXReR8vKjGqvqU1+VYiYvDFAvidxJwLC6XylUicqin7HsAn2N0CubjMMLidaC/iFSAi/IpmedFWAac5vWL7wKcTvbZEp8A/kM8R7KI7JGljFuB3in2LwR+ENsQL4d0HMtwXS+xLJCDPVk6ct7d2RHOe0qKejG59sHlgPgL8EtgRIJqjTE/RAJex6VLBViK+z/1FJc/4hRfvd7ABq+dc+LamIPr2rvTk6kI2FddFsefeNe0q1f3YDo7Qmw3xhSHEQrq0mh+A7hRnKN7FfC1DNt4AffF+izwDC467sosRbod50tY7clzdpYyVuG6olbFHMpx/BAo95y7rwJTE9T5X6DI61a7Dzhf22bVS8Rc4Mee8/uABPuvAx4QkeeBTe20Bc4/9KyIrAKuxfkf4rkNd78SOcf/zg7r4QXvOl7E5fJ+zlfvatz/7ingn3Ft3I3zzcTSC0SAv3j3ZSUuEdfH3r7R3jmNTsCi4xqGkXM8X8ocVR3XgTa+AUxU1XPbqbcXcI+qjs32XEZmmI/DMIycoy5D4x9FZDfNYi6HiPwWOIEkc2LiGIwbiWZ0EmZxGIZhGBlhPg7DMAwjI0xxGIZhGBlhisMwDMPICFMchmEYRkaY4jAMwzAy4v8Ds8NLTMdouL0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot best-fit result\n",
    "\n",
    "best_epo = best_fit.params[\"epo\"]\n",
    "best_per = best_fit.params[\"per\"]\n",
    "best_rprs = best_fit.params[\"rprs\"]\n",
    "best_ars = best_fit.params[\"ars\"]\n",
    "best_b = best_fit.params[\"b\"]\n",
    "best_u1 = best_fit.params[\"u1\"]\n",
    "best_u2 = best_fit.params[\"u2\"]\n",
    "\n",
    "phase = (x - best_epo) % best_per\n",
    "phase[phase > 0.5*best_per] -= best_per\n",
    "idx_sort = np.argsort(phase)\n",
    "\n",
    "model = transit_model(x, best_epo, best_per, best_rprs, best_ars, best_b, best_u1, best_u2)\n",
    "\n",
    "plt.plot(phase, y, \"k.\")\n",
    "plt.plot(phase[idx_sort], model[idx_sort], \"r\")\n",
    "plt.xlabel(\"Time from centre of transit (days)\")\n",
    "plt.ylabel(\"Relative flux\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
